#include<stdio.h>
#include<Arcroll.h>
// #include<Algorithm/delaunay.h>
#include<Mesh_IO/Mesh_IO.h>
#include<Mesh/libcell_iterator.h>
#include<Algorithm/marching_cubes.h>
#include<Algorithm/points_cloud_algorithm.h>
#include<tool/libcell_kd_node_iterator.h>
// #include<tool/libcell_octree_iterator.h>
#include <tool/libcell_tools_kd_node.h>
#include<Algorithm/poisson_construct.h>
#include <tool/libcell_tools_view.h>
#include<Eigen/IterativeLinearSolvers>
#include<Eigen/Sparse>
#include <tool/viewer_tools.h>
#include<vector>
#include <tool/libcell_tools_algorithm.h>
#include<libcell_macros.h>
#include <Math/LB_Quaternions.h>

using namespace Eigen;


#include <unordered_map>

typedef struct Greedy_Projection_Vert_Prop{
    double normal[3];
    double p[3];
    double theta;
    int is_legal; 
    int is_init; 
    int active;

}Greedy_Projection_Vert_Prop;
static inline void greedy_projection_vert_prop_init(Greedy_Projection_Vert_Prop* gpvp)
{
    memset(gpvp->normal,0,sizeof(double)*3);
    memset(gpvp->p,0,sizeof(double)*3);

    gpvp->is_legal=0;
    gpvp->theta=0;
    gpvp->is_init=1;
    gpvp->active=0;

}

typedef struct MY_HashKey :public std::pair<int, int> {
    MY_HashKey() = default;
    MY_HashKey(int a, int b)
    {
        this->first = a;
        this->second = b;
    }
    std::size_t operator()(const MY_HashKey& s) const noexcept
    {
        const std::size_t h1 = std::hash<int>{}(s.first);
        const std::size_t h2 = std::hash<int>{}(s.second);
        return h1 & (h2 << 1);
    }

}MY_HashKey;

static inline double compute_visual_avarage_distance(Node* n)
{
    double* box=mesh_compute_bounding_box_from_node(n);
    double min1=box[3]-box[0],min2=box[4]-box[1],min3=box[5]-box[2];
    double mid=min1>=min2?(min3>=min1?(min1):(min2>=min3?(min2):(min3))):(min3>=min2?(min2):(min1>=min3?(min1):(min3)));
    double vdis=mid/70.0;
    return vdis;
}
// static void fill_normal_coeff(Poisson_Grid*pg,Node* node_vertices)
// {
//     int i=0,j=0,k=0,id=0;
//     double* normal=NULL;
//     for(Node* nit=node_vertices;nit!=NULL;nit=(Node*)(nit->Next))
//     {
//         template_v* v=(template_v*)(nit->value);
//         i=(int)((v->point[0]+pg->range)/pg->width);
//         j=(int)((v->point[1]+pg->range)/pg->width);
//         k=(int)((v->point[2]+pg->range)/pg->width);
//         i=(i<0?0:(i>=pg->size? (pg->size-1): i ));
//         j=(j<0?0:(j>=pg->size? (pg->size-1): j )); 
//         k=(k<0?0:(k>=pg->size? (pg->size-1): k ));
//         if(pg->pgps[i][j][k]!=NULL)
//         {
//             continue;
//         }
//         Poisson_Grid_Prop*pgp=(Poisson_Grid_Prop*)malloc(sizeof(Poisson_Grid_Prop));
//         pg->pgps[i][j][k]=pgp;
//         poisson_grid_prop_init(pgp);
//         pgp->id=id;
//         pgp->pos[0]=(i+0.5)*pg->width-pg->range;
//         pgp->pos[1]=(j+0.5)*pg->width-pg->range;
//         pgp->pos[2]=(k+0.5)*pg->width-pg->range;
//         pgp->indexs[0]=i;pgp->indexs[1]=j;pgp->indexs[2]=k;
//         pg->candidate=node_overlying(pg->candidate,pgp);
        
//         normal=(double*)(v->prop);
//         pgp->normal_coeff[0]=normal[0];
//         pgp->normal_coeff[1]=normal[1];
//         pgp->normal_coeff[2]=normal[2];
        
//         id++;
//     }

//     for(Node* nit=pg->candidate;nit!=NULL;nit=(Node*)(nit->Next))
//     {
//         Poisson_Grid_Prop*pgp=(Poisson_Grid_Prop*)(nit->value);
//         int i=pgp->indexs[0],j=pgp->indexs[1],k=pgp->indexs[2];

//         for(int l=0;l<5;l++)
//         {
//             for(int m=0;m<5;m++)
//             {
//                 for(int n=0;n<5;n++)
//                 {
//                     if(i+l-2<0||i+l-2>=pg->size||j+m-2<0||j+m-2>=pg->size||k+n-2<0||k+n-2>=pg->size )
//                     {
//                         continue;
//                     }

//                     if(pg->pgps[i+l-2][j+m-2][k+n-2]==NULL)
//                     {

//                         Poisson_Grid_Prop*pgp1=(Poisson_Grid_Prop*)malloc(sizeof(Poisson_Grid_Prop));
//                         pg->pgps[i+l-2][j+m-2][k+n-2]=pgp1;
//                         poisson_grid_prop_init(pgp1);
//                         pgp1->id=id;
//                         pgp1->pos[0]=(i+l-2+0.5)*pg->width-pg->range;
//                         pgp1->pos[1]=(j+m-2+0.5)*pg->width-pg->range;
//                         pgp1->pos[2]=(k+n-2+0.5)*pg->width-pg->range;

//                         pgp1->indexs[0]=i+l-2;pgp1->indexs[1]=j+m-2;pgp1->indexs[2]=k+n-2;

//                         pg->candidate=node_overlying(pg->candidate,pgp1);

//                         // normal=(double*)(v->prop);
//                         // pgp1->normal_coeff[0]=normal[0];
//                         // pgp1->normal_coeff[1]=normal[1];
//                         // pgp1->normal_coeff[2]=normal[2];
                        
//                         id++;
//                     }    
//                 }
//             }
//         }

//     }

// }


static void test_show_triangles(Viewer_World* vw, Node* node_tri)
{
    Node* n=vw->create_something(vw,"Faces");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Faces* vf=(Viewer_Faces*)(vs->evolution);
     
    vf->Data_rows=node_size(node_tri)*3;
    vf->Data_index_rows=node_size(node_tri);
    printf("node_tr size:%d\n",vf->Data_rows);

    vf->Data=(float*)malloc(sizeof(float)*3*vf->Data_rows);
    vf->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*4*vf->Data_index_rows);
    vf->normal_rows=vf->Data_index_rows;
    vf->normal=(float*)malloc(sizeof(float)*3*vf->normal_rows);
    memset(vf->normal,0,sizeof(float)*3*vf->normal_rows);
    int i=0;
    for(Node* nit=node_tri;nit!=NULL;nit=(Node*)(nit->Next))
    {
        TRIANGLE* triangle=(TRIANGLE*)(nit->value);
        normalize(triangle->normals[0],3);
        for(int j=0;j<3;j++)
        {
            for(int k=0;k<3;k++)
            {
                vf->Data[i*9+j*3+k]=triangle->p[j][k];  
            }
            vf->normal[i*3+j]=triangle->normals[0][j];
        }
        vf->Data_index[i*4+0]=3; vf->Data_index[i*4+1]=i*3+0;vf->Data_index[i*4+2]=i*3+1;vf->Data_index[i*4+3]=i*3+2;
        i++;
    } 
    vf->color_rows=vf->Data_index_rows;
    vf->set_color(vf,0.9,0.6,0.7,1.0);
    free_node(n); 
}
static void fill_sparse_matrix(Eigen:: SparseMatrix<double>& H,std::unordered_map<MY_HashKey,double,MY_HashKey> & unmap)
{
    std::vector<Eigen::Triplet<double>> tripletList;
    for(auto it =unmap.begin();it!=unmap.end();it++)
    {
        tripletList.push_back(Eigen::Triplet<double>(it->first.first,it->first.second,it->second));
    }
    H.setFromTriplets(tripletList.begin(),tripletList.end());

}
// static void  poisson_test_sparse_matrix(Poisson_Grid*pg )
// {
//     int len=node_size(pg->candidate);

//     SparseMatrix<double,RowMajor> A(len,len);
//     A.reserve(VectorXi::Constant(len,140));
//     printf("begin sparse matrix %d \n",len);
//     //std::unordered_map<MY_HashKey,double,MY_HashKey> unmap;

//     int id=0,id1=0;
//     int sum=0;double temp=0;

//     for(Node* nit= pg->candidate;nit!=NULL;nit=(Node*)(nit->Next))
//     {
//         Poisson_Grid_Prop*pgp=(Poisson_Grid_Prop*)(nit->value);
//         id=pgp->id;
//         sum++;
//         int i=pgp->indexs[0],j=pgp->indexs[1],k=pgp->indexs[2];

//         for(int l=0;l<5;l++ )
//         {
//             for(int m=0;m<5;m++)
//             {
//                 for(int n=0;n<5;n++)
//                 {
//                     if(i+l-2<0||i+l-2>=pg->size||j+m-2<0||j+m-2>=pg->size||k+n-2<0||k+n-2>=pg->size)
//                     {
//                         continue;
//                     }
//                     if(pg->pgps[i+l-2][j+m-2][k+n-2]==NULL)
//                     {
//                         continue;
//                     }
//                     id1=pg->pgps[i+l-2][j+m-2][k+n-2]->id;
//                     temp=0;
//                     temp+=(fx_dot_DDfx[l]*fx_dot_fx[m]*fx_dot_fx[n]); 
//                     temp+=(fx_dot_fx[l]*fx_dot_DDfx[m]*fx_dot_fx[n]); 
//                     temp+=(fx_dot_fx[l]*fx_dot_fx[m]*fx_dot_DDfx[n]);
//                     temp/=(pg->width*pg->width);
//                     A.coeffRef(id,id1)=temp;

//                 }
//             }
//         }
//     }      
    
//     printf("beigin compress\n");
//     //fill_sparse_matrix(A,unmap);
//     A.makeCompressed();

//     ConjugateGradient<SparseMatrix<double>,Eigen::Upper> solver;
//     //solver.setMaxIterations(50);
//     //Eigen::BiCGSTAB<Eigen::SparseMatrix<double>, Eigen::IncompleteLUT<double>> solver;
//     //Eigen::SimplicialLDLT<Eigen::SparseMatrix<double>>solver;
//     double*bi= poisson_compute_b(pg );

//     Eigen::VectorXd b(len);
//     for(int i=0;i<len;i++)
//     {
//         b.coeffRef(i)=bi[i];
//     }   

//     printf("end b\n"); 

//     free(bi);
      
//     VectorXd x= solver.compute(A).solve(b); 
//     pg->solve=(double*)malloc(sizeof(double)*len); 
//     for(int i=0;i<len;i++)
//     {
//         pg->solve[i]=x.coeffRef(i);
//     }
//     printf(" 1sum:%d\n",sum);
// }
// void  test_show_grid_normal(Viewer_World*vw,Poisson_Grid*pg)
// { 
//     Node*n2=vw->create_something(vw,"Edges");
//     Viewer_Something* vs=(Viewer_Something*)(n2->value);
//     Viewer_Edges* ve=(Viewer_Edges*)(vs->evolution);
//     int sum=0;

   
//         if(pg->pgps[i][j][k].normal_coeff[0]!=0&&pg->pgps[i][j][k].normal_coeff[1]!=0&&pg->pgps[i][j][k].normal_coeff[2]!=0) 
//         {
//             sum++;
//         }
//     ve->Data_rows=sum*2;
//     ve->Data_index_rows=sum;
//     ve->color_rows=ve->Data_index_rows;
//     ve->Data=(float *)malloc(sizeof(float)*3*ve->Data_rows);    
//     ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows);
//     ve->set_color(ve,0,1.0,0.0,1.0);  
//     sum=0;
//     //printf("data rows:%d index rows:%d\n",ve->Data_rows,ve->Data_index_rows);
//     for(int i=0;i<pg->size;i++)
//     {
//         for(int j=0;j<pg->size;j++)
//         {
//             for(int k=0;k<pg->size;k++)
//             {
//                 if(pg->pgps[i][j][k].normal_coeff[0]!=0&&pg->pgps[i][j][k].normal_coeff[1]!=0&&pg->pgps[i][j][k].normal_coeff[2]!=0) 
//                 {
//                     ve->Data[(sum*2)*3+0]=pg->pgps[i][j][k].pos[0];
//                     ve->Data[(sum*2)*3+1]=pg->pgps[i][j][k].pos[1];
//                     ve->Data[(sum*2)*3+2]=pg->pgps[i][j][k].pos[2];
                    
//                     ve->Data[(sum*2+1)*3+0]= ve->Data[(sum*2)*3+0] +pg->pgps[i][j][k].normal_coeff[0]*0.1;
//                     ve->Data[(sum*2+1)*3+1]= ve->Data[(sum*2)*3+1] +pg->pgps[i][j][k].normal_coeff[1]*0.1;
//                     ve->Data[(sum*2+1)*3+2]= ve->Data[(sum*2)*3+2] +pg->pgps[i][j][k].normal_coeff[2]*0.1;
                    
//                     //printf("%lf %lf %lf %lf %lf %lf\n",pg->pgps[i][j][k].pos[0],pg->pgps[i][j][k].pos[1],pg->pgps[i][j][k].pos[2],pg->pgps[i][j][k].normal_coeff[0],pg->pgps[i][j][k].normal_coeff[1],pg->pgps[i][j][k].normal_coeff[2]);
//                     ve->Data_index[sum*2]=sum*2;
//                     ve->Data_index[sum*2+1]=sum*2+1; 
//                     sum++;
//                     //vp->Data[sum*3+0]=pg->pgps ; 
//                 }
//             }
//         }
//     }  
//     free_node(n2);
// }


// static inline int is_illegal_tri(TRIANGLE* t)
// {
//     if(is_same_point(t->p[0],t->p[1])||is_same_point(t->p[2],t->p[1])||is_same_point(t->p[0],t->p[2]))
//     {
//         return 0;
//     }
//     return 1;
// }



Mesh* create_mesh_from_marching_cubesn(Node****tris,Node* tris1,int len_grids)
{
    Mesh* re=(Mesh*)malloc(sizeof(Mesh));
    Mesh_init(re);
    re->dimension=2;re->simplex=1;re->manifold_require=1;
    for(Node*nit =tris1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        TRIANGLE*t=(TRIANGLE*)(nit->value);
        if(t->prop==NULL)
        {
            template_v*vs[3]={NULL};
            vs[0]=re->create_vertexv(re,t->p[0],3);
            vs[1]=re->create_vertexv(re,t->p[1],3);
            vs[2]=re->create_vertexv(re,t->p[2],3);
            template_c*c= re->create_cellv(re,vs,3);
            t->prop=c;
        }
    }
    return re;

}


static Viewer_Something* test_show_mesh_cells(Viewer_World* vw,Mesh*m)
{
    Node* n=vw->create_something(vw,"Faces");
    Viewer_Something*vs=(Viewer_Something*)(n->value);
    Viewer_Faces* vf=(Viewer_Faces*)(vs->evolution);
     
    vf->Data_rows=m->num_v(m);
    vf->Data_index_rows=m->num_c(m);
    get_data_from_2dim_cell(m,&(vf->Data),&(vf->Data_index));
    

    
    vf->color_rows=vf->Data_index_rows;
    vf->set_color(vf,0.4,0.6,0.7,1.0);

    vf->normal_rows=vf->Data_index_rows;
    vf->compute_normal(vf);
    
    free_node(n); 
    return vs;

}

// void test_poisson()
// {
//     Mesh mesh;
//     Mesh_init(&mesh);
//     _ReadOff_(&mesh,"rabbitCut.off",3);
//     //_ReadPly_(&mesh,"bunny.points.ply");
//     mesh_adjust_mesh_to_unit_hexahedron(&mesh); 

//     template_v ** values=(template_v**)malloc(sizeof(template_v*)*mesh.num_v(&mesh));
//     int i=0;
//     Node* node_vertices=NULL;
//     for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++ )
//     {
//         values[i]=quote(vit);
//         node_vertices=node_overlying(node_vertices,quote(vit));
//         i++;
//     }
//     KD_Node*kdtree1= create_kd_tree(values,i,0); 
//     points_cloud_compute_normal(node_vertices,kdtree1,8);
//     double r=compute_visual_avarage_distance(node_vertices);
//     adjust_normal_of_node_vertices(node_vertices,r,kdtree1);
    
//     Poisson_Grid*pg=(Poisson_Grid*)malloc(sizeof(Poisson_Grid));
//     poisson_grid_init(pg);

//     poisson_grid_init_with_depth(pg,7);

//     fill_normal_coeff(pg,node_vertices);

//     poisson_test_sparse_matrix(pg);
    
//     double iso=compute_average_poisson_iso(pg);
//     printf("begin hrerere\n");
//     // Node* node_tris=marching_cubes_generate_triangles(0.2,100,my_fun2,NULL);
//     // Node* node_tris=marching_cubes_generate_triangles(iso,100,end_poisson_space_fun,pg);
//     Node*node_tris1=NULL;
//     Node****node_tris=marching_cubes_generate_trianglesn(iso,100,end_poisson_space_fun,pg,&node_tris1);

//     Mesh* m1= create_mesh_from_marching_cubes(node_tris,node_tris1,100);
//     // template_f *f=m1->get_facep(m1,1522);
//     // if(f->halffaces[0]->cell!=NULL&&f->halffaces[1]->cell!=NULL)
//     // {

//     //     printf("cuowuwuuwuwuwuwu********\n");
//     // }


//     template_c* c1=m1->get_cellp(m1,40583);
//     template_c* c2=m1->get_cellp(m1,40584);

//     printf("dd:%d %d %d\n",c1->vertices[0]->id,c1->vertices[1]->id,c1->vertices[2]->id);
    
//     printf("%lf %lf %lf\n",c1->vertices[0]->point[0],c1->vertices[0]->point[1],c1->vertices[0]->point[2] );
//     printf("%lf %lf %lf\n",c1->vertices[1]->point[0],c1->vertices[1]->point[1],c1->vertices[1]->point[2] );
//     printf("%lf %lf %lf\n",c1->vertices[2]->point[0],c1->vertices[2]->point[1],c1->vertices[2]->point[2] );


//     printf("dd:%d %d %d\n",c2->vertices[0]->id,c2->vertices[1]->id,c2->vertices[2]->id);

//     printf("%lf %lf %lf\n",c2->vertices[0]->point[0],c2->vertices[0]->point[1],c2->vertices[0]->point[2] );
//     printf("%lf %lf %lf\n",c2->vertices[1]->point[0],c2->vertices[1]->point[1],c2->vertices[1]->point[2] );
//     printf("%lf %lf %lf\n",c2->vertices[2]->point[0],c2->vertices[2]->point[1],c2->vertices[2]->point[2] );


//     Mesh*m2= create_mesh_from_marching_cubesn(node_tris,node_tris1,100);

//     //create_mesh_from_marching_cubesn(m1,node_tris,node_tris1,100);
//     //_WriteOff_(m1,"rabit.off");
//     Viewer_World_Manager vwm;
//     viewer_world_manager_init(&vwm);
//     Viewer_World *vw=vwm.create_world(&vwm,NULL);

//     add_default_somethings(vw);  
//    // test_show_vertices(vw,vertices1);
//     //test_show_grid_points(vw,pg);
//     //test_show_grid_normal(vw,pg); 
//     //test_show_mesh_cells(vw,m1);
//     // Viewer_Something* vs= test_show_mesh_cells(vw,m2);
//     // Viewer_Faces* vf=(Viewer_Faces*)(vs->evolution);
//     // //float color[4]={1,0,0,1};
//     // vf->set_color(vf,1,0,0,1);
//     //test_show_triangles(vw, Node* node_tri)
//     test_show_mesh_cells(vw,m1);
//     Node* n=vw->create_something(vw,"Points");
//     Viewer_Something*vs=(Viewer_Something*)(n->value);
//     Viewer_Points* vp=(Viewer_Points*)(vs->evolution);
     
//     vp->Data_rows=3;
//     vp->Data=(float*)malloc(sizeof(float)*3*vp->Data_rows);
//     viewer_tools_fill_data("%f",vp->Data_rows*3,&(vp->Data),
//         -0.44,0.6,-0.06,
//         -0.42,0.6,-0.04,
//         -0.44,0.6,-0.06); 
    
//     vp->set_color(vp,1.0,0.0,0.0,1.0);

//     free_node(n); 

//     //test_show_triangles(vw,node_tris1);
 
//     Viewer_Opengl_Interpreter voi;
//     viewer_opengl_interpreter_initn(&voi,&vwm);
//     voi.interpreter(&voi);
//     viewer_opengl_interpreter_free(&voi);
//     viewer_world_manager_free(&vwm);
//     poisson_grid_free(pg);
 
//     free_kdnode(kdtree1);
//     free_node(node_vertices);
//     free(values);


//     Mesh_free(&mesh);
// }

static inline double* compute_triangle_normal(double* q1,double* q2,double* q3)
{
    double p1[3]={q2[0]-q1[0],q2[1]-q1[1],
        q2[2]-q1[2]};
    double p2[3]={q3[0]-q1[0],q3[1]-q1[1],
        q3[2]-q1[2]};
    double*re=out_product(p1,p2); 
    normalize(re,3);
    return re;
}


static double* compute_cell_normal(template_c*c)
{
    if(c==NULL)
    {
        return NULL;
    }
    template_hf* hf=(template_hf*)(c->halffaces->value);
    template_v* v=NULL;
    for(int i=0;i<c->vertices_size;i++)
    {
        if(c->vertices[i]!=hf->vertices[0]&&c->vertices[i]!=hf->vertices[1])
        {
            v=c->vertices[i];
            break;
        }
    }
    if(v==NULL)
    {return NULL;}
    
    double p1[3]={hf->vertices[0]->point[0]-v->point[0],hf->vertices[0]->point[1]-v->point[1],
        hf->vertices[0]->point[2]-v->point[2]};

    double p2[3]={hf->vertices[1]->point[0]-v->point[0],hf->vertices[1]->point[1]-v->point[1],
        hf->vertices[1]->point[2]-v->point[2]};

    double*re=out_product(p1,p2); 
    normalize(re,3);
    return re;
}

static int is_legal_triangle_angle(double *p1,double* p2,double* p3 );

//***************
//大于0是在外心圆内部
//小于0是在外心圆外部
//****************
static inline double is_concyclic_four_vertices(double* p1,double* p2,double* p3,double* p4)
{
    LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
    lb_matrix_init_double(mat,4,4); 
    double* data=(double*)(mat->data);
    data[0*4+0]=p1[0]; data[0*4+1]=p1[1];data[0*4+2]=p1[0]*p1[0]+p1[1]*p1[1];data[0*4+3]=1.0;
    data[1*4+0]=p2[0]; data[1*4+1]=p2[1];data[1*4+2]=p2[0]*p2[0]+p2[1]*p2[1];data[1*4+3]=1.0;
    data[2*4+0]=p3[0]; data[2*4+1]=p3[1];data[2*4+2]=p3[0]*p3[0]+p3[1]*p3[1];data[2*4+3]=1.0;
    data[3*4+0]=p4[0]; data[3*4+1]=p4[1];data[3*4+2]=p4[0]*p4[0]+p4[1]*p4[1];data[3*4+3]=1.0;
    double *re=((double*)(mat->det(mat)));
    double re1=*re;
    lb_matrix_free(mat);
    free(re);
    return re1;
}

static inline double compute_triangle_area(double* p1,double* p2,double* p3)
{
    double q1[3]={p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]};
    
    double q2[3]={p3[0]-p1[0],p3[1]-p1[1],p3[2]-p1[2]};

    double q[3]={q1[1]*q2[2]-q1[2]*q2[1],-q1[0]*q2[2]+q1[2]*q2[0],q1[0]*q2[1]-q1[1]*q2[0]};
    double re=SAFE_SQRT(q[0]*q[0]+q[1]*q[1]+q[2]*q[2])/2.0;

    return re;
}

static inline double compute_score_triangle(template_v* v1,template_v* v2,template_v* v3 )
{
    double normal[3]={0};
    Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v1->prop);
    normal[0]+=1.0*gpvp->normal[0]; normal[1]+=1.0*gpvp->normal[1];normal[2]+=1.0*gpvp->normal[2];
    gpvp=(Greedy_Projection_Vert_Prop*)(v2->prop); 
    normal[0]+=gpvp->normal[0]; normal[1]+=gpvp->normal[1];normal[2]+=gpvp->normal[2];
    gpvp=(Greedy_Projection_Vert_Prop*)(v3->prop); 
    normal[0]+=gpvp->normal[0]; normal[1]+=gpvp->normal[1];normal[2]+=gpvp->normal[2]; 
    normalize(normal,3);

    double* n1=compute_triangle_normal(v1->point,v2->point,v3->point );

    //double re=n1[0]*gpvp->normal[0]+n1[1]*gpvp->normal[1]+n1[2]*gpvp->normal[2]+1.0 ;
    double re=n1[0]*normal[0]+n1[1]*normal[1]+n1[2]*normal[2];

    SAFE_FREE(n1);

    return re;
}


static inline double compute_weight_score_triangle(template_v* v1,template_v* v2,template_v* v3)
{
    double score=compute_score_triangle(v1,v2,v3);
    double weight=compute_triangle_area(v1->point,v2->point,v3->point);

    return score*weight;

}

static inline double compute_weight_score_2triangles(template_v* v1,template_v* v2,template_v* v3,
    template_v* v4,template_v* v5,template_v* v6 )
{
    double score1=compute_score_triangle(v1,v2,v3);
    double weight1=compute_triangle_area(v1->point,v2->point,v3->point);
    double score2=compute_score_triangle(v4,v5,v6);
    double weight2=compute_triangle_area(v4->point,v5->point,v6->point);
 
    return (score1*weight1+score2*weight2)/(weight1+weight2);
}
// static inline int four_vertices_is_same_dir(template_v** vs)
// {
//     for(int i=0;i<3;i++)
//     {
//         Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(vs[i]->prop);
//         Greedy_Projection_Vert_Prop*gpvp1=(Greedy_Projection_Vert_Prop*)(vs[i+1]->prop);
//         if(gpvp->normal[0]*gpvp1->normal[0]+gpvp->normal[0]*gpvp1->normal[0]+
//             gpvp->normal[0]*gpvp1->normal[0]<0.9 )
//         {
//             return 0;
//         }
//     }
//     return 1;
// }

static inline  template_f* swap_edge(Mesh* m, template_f* f,template_v** vs)
{
    printf("begin swap\n");
    template_c* c1=f->halffaces[0]->cell;
    template_c* c2=f->halffaces[1]->cell; 
    m->delete_cell(m,*c1,true);
    m->delete_cell(m,*c2,true);
    
    template_v* vs1[3]={vs[3],vs[0],vs[1]};
    m->create_cellv(m,vs1,3);
    vs1[0]=vs[1];vs1[1]=vs[2];vs1[2]=vs[3];
    m->create_cellv(m,vs1,3);

    template_v* vs2[2]={vs[1],vs[3]};
    template_f* re= m->get_facev(m,vs2,2);
     
    return re; 
}

static template_f* touying_tiaozheng_edge(Mesh* m, template_f* f)
{
    printf("once tiaozheng\n");  
    printf("once tiaozhengherer\n");

    template_v* vs[4]={0};

    template_c* c1=f->halffaces[0]->cell;
    template_c* c2=f->halffaces[1]->cell;
 
    vs[0]=f->halffaces[0]->vertices[0]; vs[2]= f->halffaces[0]->vertices[1];
    printf("once tiaozhengherer1 %d %d\n",c1->id,c2->id);
    printf("once tiaozhengherer2 %d %d %d %d %d %d\n",c1->vertices[0]->id,c1->vertices[1]->id,c1->vertices[2]->id,
        c2->vertices[0]->id,c2->vertices[1]->id,c2->vertices[2]->id);

    for(int i=0;i<3;i++)
    {
        template_v* v=c1->vertices[i];

        if(v!=vs[0]&&v!=vs[2])
        {
            vs[3]=v;
            break;
        }
    }  
    for(int i=0;i<3;i++)
    {
        template_v* v=c2->vertices[i];
        if(v!=vs[0]&&v!=vs[2])
        {
            vs[1]=v;
            break;
        }
    } 
    printf("once tiaozhengherer2 %d %d %d %d %d %d\n",c1->vertices[0]->id,c1->vertices[1]->id,c1->vertices[2]->id,
        c2->vertices[0]->id,c2->vertices[1]->id,c2->vertices[2]->id);

    Greedy_Projection_Vert_Prop*gpvps[4]={0};
    gpvps[0]=(Greedy_Projection_Vert_Prop*)(vs[0]->prop);
    gpvps[1]=(Greedy_Projection_Vert_Prop*)(vs[1]->prop);
    gpvps[2]=(Greedy_Projection_Vert_Prop*)(vs[2]->prop);
    gpvps[3]=(Greedy_Projection_Vert_Prop*)(vs[3]->prop);
    printf("is is_concyclic_four_vertices:%lf\n",
        is_concyclic_four_vertices(gpvps[0]->p,gpvps[1]->p,gpvps[2]->p,gpvps[3]->p));

    printf("%d %d %d %d\n",vs[0]->id,vs[1]->id,vs[2]->id,vs[3]->id);

    if(is_concyclic_four_vertices(gpvps[0]->p,gpvps[1]->p,gpvps[2]->p,gpvps[3]->p)<=0)
    {
        return NULL;
    }
    // double normal[3]={0};
    // normal[0]+=gpvps[0]->normal[0]; normal[1]+=gpvps[0]->normal[1];normal[2]+=gpvps[0]->normal[2];
    // normal[0]+=gpvps[1]->normal[0]; normal[1]+=gpvps[1]->normal[1];normal[2]+=gpvps[1]->normal[2];
    // normal[0]+=gpvps[2]->normal[0]; normal[1]+=gpvps[2]->normal[1];normal[2]+=gpvps[2]->normal[2];
    // normal[0]+=gpvps[3]->normal[0]; normal[1]+=gpvps[3]->normal[1];normal[2]+=gpvps[3]->normal[2];
 
    // normalize(normal,3);

    // double* normal1=compute_cell_normal(c1);

    // double* normal2=compute_cell_normal(c2);
    //compute_weight_score_cell
    // double score1=compute_score_cell(c1), score2=compute_score_cell(c2);
    // double score1=compute_weight_score_triangle(vs[1],vs[2],vs[0])+
    //             compute_weight_score_triangle(vs[3],vs[0],vs[2]);

    // double score2=compute_weight_score_triangle(vs[0],vs[1],vs[3])+
    //             compute_weight_score_triangle(vs[2],vs[3],vs[1]);

    double score1=compute_weight_score_2triangles(vs[1],vs[2],vs[0],
        vs[3],vs[0],vs[2]);

    double score2=compute_weight_score_2triangles(vs[0],vs[1],vs[3],
        vs[2],vs[3],vs[1]);
    
    printf("score1:%lf score2:%lf\n",score1,score2);
    if(score2>=0.97*score1)
    {
        return swap_edge(m,f,vs);
    } 
   
    // if(!is_legal_triangle_angle(vs[0]->point,vs[1]->point,vs[2]->point)
    //     ||!is_legal_triangle_angle(vs[3]->point,vs[0]->point,vs[2]->point) )
    // { 
    //     if(is_legal_triangle_angle(vs[3]->point,vs[0]->point,vs[1]->point)
    //         &&is_legal_triangle_angle(vs[1]->point,vs[2]->point,vs[3]->point))
    //     {
    //         //compute_score_triangle
    //         double score1=compute_weight_score_triangle(vs[0],vs[1],vs[2])+
    //             compute_weight_score_triangle(vs[3],vs[0],vs[2]);
    //         double score2=compute_weight_score_triangle(vs[3],vs[0],vs[1])+
    //             compute_weight_score_triangle(vs[1],vs[2],vs[3]);
    //         if(score2>=score1)
    //         {
    //             swap_edge(m,f,vs);
    //             return 1;
    //         }
    //     } 
    // }

 
    return NULL;
}



typedef struct Cell_Prop{
    double* normal;
    Node* node;

}Cell_Prop;
static inline void cell_prop_init(Cell_Prop*cp)
{
    cp->node=NULL;
    cp->normal=NULL;
}


static void my_points_cloud_compute_normalr(Node* n,KD_Node*kdtree,double r)
{
    int i=0;
    double* normal=NULL;
    for(Node* nit=n;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        Node* n1=kd_tree_find_nearest_sphere(v,r,kdtree);
        double**axis=compute_pca_result(n1);
        template_v* v1=(template_v*)(v->prop);
        template_c*c=(template_c*)(v1->prop);
        Cell_Prop* cp=(Cell_Prop*)(c->prop);


        normal=(double*)malloc(sizeof(double)*3);
        v->prop=normal;
        if(axis==NULL)
        {
            normal[0]=0;normal[1]=0;normal[2]=1.0;
        }
        else
        {
            normal[0]=axis[0][0];normal[1]=axis[1][0];normal[2]=axis[2][0];
            for(int j=0;j<v->point_size;j++)
            {
                free(axis[j]);
            }
            free(axis);
        }
        if(inner_product(cp->normal,normal,3)<0)
        {
            normal[0]=-normal[0];normal[1]=-normal[1];normal[2]=-normal[2];
        }
        free_node(n1);
        i++;
    } 


}

static double* marching_cubes_compute_normal(double* p,double step,double ( * fun)(double* p,void* tool),void* tool1 )
{
    double* re=(double*)malloc(sizeof(double)*3);
    double temp_p[3]={0};
    double value=fun(p,tool1);
    for(int i=0;i<3;i++)
    {
        memmove(temp_p,p,sizeof(double)*3);
        temp_p[i]+=step;
        double value1=fun(temp_p,tool1);
        re[i]=(value1-value)/step;  
    }
    return re;
}

static inline template_v* compute_hou_point_from_vertex(Mesh* m,template_v*v,double hou ,Mesh*m1)
{
    double p[3]={0};
    double normal[3]={0};
    for(auto vcit=m->vc_begin(m,*v);vcit!=m->vc_end(m,*v);vcit++)
    {
        Cell_Prop* cp=(Cell_Prop*)(quote(vcit)->prop);
        if(cp!=NULL)
        {
            normal[0]+=cp->normal[0];
            normal[1]+=cp->normal[1];
            normal[2]+=cp->normal[2];
        }
    }
    template_v *v1=NULL,*v2=v,*v3=NULL;

    for(auto vfit=m->vf_begin(m,*v);vfit!=m->vf_end(m,*v);vfit++)
    {
        if(m->face_is_boundary(m,*vfit))
        {
            template_hf* hf=NULL;
            if(quote(vfit)->halffaces[0]->cell!=NULL)
            {
                hf=quote(vfit)->halffaces[0];
            }
            else
            {
                hf=quote(vfit)->halffaces[1];
            }
            if(hf->vertices[0]==v2)
            {
                v3=hf->vertices[1];
            }
            else if(hf->vertices[1]==v2)
            {
                v1=hf->vertices[0];
            }
        } 
    }

    if(normalize(normal,3))
    {
        double p1[3]={v1->point[0]-v2->point[0],v1->point[1]-v2->point[1],v1->point[2]-v2->point[2]};
        double p2[3]={v3->point[0]-v2->point[0],v3->point[1]-v2->point[1],v3->point[2]-v2->point[2]};
        double inn1=inner_product(p1,normal,3);
        double inn2=inner_product(p2,normal,3);
        p1[0]-=normal[0]*inn1;p1[1]-=normal[1]*inn1;p1[2]-=normal[2]*inn1;
        p2[0]-=normal[0]*inn2;p2[1]-=normal[1]*inn2;p2[2]-=normal[2]*inn2;
        normalize(p1,3); normalize(p2,3);

        double p4[3]={(p1[0]+p2[0])/2.0,(p1[1]+p2[1])/2.0,(p1[2]+p2[2])/2.0};

        double* p3=out_product(p1,p2);

        if(inner_product(normal,p3,3)>0)
        {
            p4[0]=-p4[0];p4[1]=-p4[1];p4[2]=-p4[2];
        }
        normalize(p4,3);

        SAFE_FREE(p3);
        p[0]=v2->point[0]+hou*normal[0];        
        p[1]=v2->point[1]+hou*normal[1];
        p[2]=v2->point[2]+hou*normal[2]; 
        template_v* re=m1->create_vertexv(m1,p,3);
        double* normal1=(double*)malloc(sizeof(double)*3);
        memmove(normal1,normal,sizeof(double)*3);
        re->prop=normal1;

        return re;
    }

    return NULL;
}

static inline double my_atan2(double a,double b)
{
    double re= atan2(a,b);
    a<0?(re=re+2*M_PI):(re=re);
    return re;
}


static int is_vertex_in_cell(template_v*v,template_c* c)
{
    if(v==c->vertices[0]||v==c->vertices[1]||v==c->vertices[2])
    {
        return 0;
    }
    Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
    double ps[3]={0};
    int i=0;
    for(Node* nit=c->halffaces;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf* hf=(template_hf*)(nit->value);
 
        Greedy_Projection_Vert_Prop*gpvp1=(Greedy_Projection_Vert_Prop*)(hf->vertices[0]->prop);
        Greedy_Projection_Vert_Prop*gpvp2=(Greedy_Projection_Vert_Prop*)(hf->vertices[1]->prop);
        double p1[3]={gpvp1->p[0]-gpvp->p[0],gpvp1->p[1]-gpvp->p[1],0 };
        double p2[3]={gpvp2->p[0]-gpvp->p[0],gpvp2->p[1]-gpvp->p[1],0 };
        ps[i]=p1[0]*p2[1]-p1[1]*p2[0];
        //ps[i]=out_product(p1,p2);
        i++;
    }
    int re=0;
    if(ps[0]*ps[1]>0&&ps[0]*ps[2]>0)
    {
        re=1;
    } 
    //SAFE_FREE(ps[0]);SAFE_FREE(ps[1]);SAFE_FREE(ps[2]);
    return re;

}
static template_c* vertex_is_in_which_cell(Node* node_c,template_v* v)
{
    for(Node* nit=node_c;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        if(is_vertex_in_cell(v, c))
        {
            return c;
        }
    }
    return NULL;
}

static int is_legal_triangle_angle(double *p1,double* p2,double* p3 )
{
    double q1[3]={p2[0]-p1[0],p2[1]-p1[1],p2[2]-p1[2]};
    double q2[3]={p3[0]-p2[0],p3[1]-p2[1],p3[2]-p2[2]};
    double q3[3]={p1[0]-p3[0],p1[1]-p3[1],p1[2]-p3[2]};
    
    double *ps[3]={q2,q3,q1};
    double dis[3]={
        SAFE_SQRT(ps[0][0]*ps[0][0]+ps[0][1]*ps[0][1]+ps[0][2]*ps[0][2] ),
        SAFE_SQRT(ps[1][0]*ps[1][0]+ps[1][1]*ps[1][1]+ps[1][2]*ps[1][2]),
        SAFE_SQRT(ps[2][0]*ps[2][0]+ps[2][1]*ps[2][1]+ps[2][2]*ps[2][2] )
    };

    int mark_min=0,mark_max=0;
    
    dis[1]>dis[mark_max]?(mark_max=1):(mark_max=mark_max);
    dis[1]<dis[mark_min]?(mark_min=1):(mark_min=mark_min);


    dis[2]>dis[mark_max]?(mark_max=2):(mark_max=mark_max);
    dis[2]<dis[mark_min]?(mark_min=2):(mark_min=mark_min); 
 
    double cos1=(dis[(mark_max+2)%3]*dis[(mark_max+2)%3]+dis[(mark_max+1)%3]*dis[(mark_max+1)%3]-
        dis[mark_max]*dis[mark_max])/(2*dis[(mark_max+1)%3]*dis[(mark_max+2)%3]);

    double cos2=(dis[(mark_min+2)%3]*dis[(mark_min+2)%3]+dis[(mark_min+1)%3]*dis[(mark_min+1)%3]-
        dis[mark_min]*dis[mark_min])/(2*dis[(mark_min+1)%3]*dis[(mark_min+2)%3]);

    printf("cos1:%lf cos2 :%lf\n",cos1,cos2);

    if(cos1>-0.95&&cos2<0.95)
    {
        return 1;
    }
    return 0;

}

static template_c* vertex_is_in_which_cell_legal(Node* node_c,template_v* v)
{
    template_c* re=NULL;
    for(Node* nit=node_c;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        if(is_vertex_in_cell(v, c))
        {
            re=c;
            break;
        }
    }
    
    return re;
}

static template_c* vertex_is_in_which_cell1(Int_RB_Tree* tree_c,template_v* v)
{
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        if(is_vertex_in_cell(v, c))
        {
            return c;
        }
    }
    return NULL;
}

static inline void touying_one_vertex(double*matd,template_v* v,template_v* v1)
{
    Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v1->prop);
    gpvp->is_init=0;
    double p[3]={v1->point[0]-v->point[0],v1->point[1]-v->point[1],v1->point[2]-v->point[2]};
    gpvp->p[0]=(matd[0*3+0]*p[0])+(matd[0*3+1]*p[1])+(matd[0*3+2]*p[2]);
    gpvp->p[1]=(matd[1*3+0]*p[0])+(matd[1*3+1]*p[1])+(matd[1*3+2]*p[2]);
    gpvp->p[2]=(matd[2*3+0]*p[0])+(matd[2*3+1]*p[1])+(matd[2*3+2]*p[2]);

}
// 0 点落在区域外
// 1 点落在区域内
static inline int is_point_in_line_segment(double*l,double*p)
{
    double data[4]={0};
    data[0*2+0]=l[0];data[0*2+1]=l[1];
    data[1*2+0]=-l[1];data[1*2+1]=l[0];

    double temp_p[2]={p[0]-l[2],p[1]-l[3]};
    double temp_p1[2]={data[0]*temp_p[0]+data[1]*temp_p[1],data[2]*temp_p[0]+data[3]*temp_p[1] };
    return (temp_p1[1]>0);
}

//0 没有交集
//1 有交集,有交叉
//2 有交集，完全落在区域内

static int  cut_line_segment_from_line(double *x,double*l,double*xn)
{
    memmove(xn,x,sizeof(double)*4);

    int re=0;
    LB_Matrix* mat=(LB_Matrix*)malloc(sizeof(LB_Matrix));
    lb_matrix_init_double(mat,2,2);
    double *data=(double*)(mat->data);
    data[0*2+0]=l[0];data[0*2+1]=-l[1];
    data[1*2+0]=l[1];data[1*2+1]=l[0];

    LB_Matrix *mati=lb_matrix_get_transpose_double(mat);

    LB_Matrix* temp_mat=(LB_Matrix*)malloc(sizeof(LB_Matrix)); 
    lb_matrix_init_double(temp_mat,2,2);
    double *temp_data=(double*)(temp_mat->data);
    temp_data[0]=x[0]-l[2];temp_data[1]=x[2]-l[2];
    temp_data[2]=x[1]-l[3];temp_data[3]=x[3]-l[3];
    
    temp_mat->mult1(mati,temp_mat);
    temp_data=(double*)(temp_mat->data);

    if(temp_data[2]<=0&&temp_data[3]<=0)
    {
        re=0; 
    }
    else if(temp_data[2]>=0&&temp_data[3]>=0)
    {
        re=2;
    }
    else if(temp_data[2]>0)
    {
        re=1;
        double lambda=temp_data[3]/(temp_data[3]-temp_data[2]);
        double p[2]={lambda*temp_data[0]+(1-lambda)*temp_data[1], lambda*temp_data[2]+(1-lambda)*temp_data[3]};
        xn[2]=data[0]*p[0]+data[1]*p[1];
        xn[3]=data[2]*p[0]+data[3]*p[1];
        xn[2]+=l[2];xn[3]+=l[3];
    }
    else if(temp_data[3]>0)
    {
        re=1;
        double lambda=temp_data[3]/(temp_data[3]-temp_data[2]);
        double p[2]={lambda*temp_data[0]+(1-lambda)*temp_data[1], lambda*temp_data[2]+(1-lambda)*temp_data[3]};
        xn[0]=data[0]*p[0]+data[1]*p[1];
        xn[1]=data[2]*p[0]+data[3]*p[1];
        xn[0]+=l[2];xn[1]+=l[3]; 
    }

    lb_matrix_free(temp_mat);
    lb_matrix_free(mati);
    lb_matrix_free(mat);
    return re;
}

static inline int my_panduan2_point_is_in(double*l[2],double* p )
{
    return (is_point_in_line_segment(l[0],p)==0||is_point_in_line_segment(l[1],p)==0)?0:1;
}

//0 无交集
//1 有交集,有交叉
//2 有交集，完全落在区域内

static inline int my_panduan3_is_in(double l[3][4],double* x)
{
    int re=0;
    double xn[4]={0};
    int temp_re=0;
    temp_re=cut_line_segment_from_line(x,l[0],xn);
    //printf("panduan1 :%lf %lf %lf %lf\n",xn[0],xn[1],xn[2],xn[3]);

    if(temp_re==0){  memmove(x,xn,sizeof(double)*4); return 0;}
    re=temp_re;
    double xn1[4]={0}; 
    temp_re=cut_line_segment_from_line(xn,l[1],xn1);
    //printf("panduan2 :%lf %lf %lf %lf\n",xn1[0],xn1[1],xn1[2],xn1[3]);

    if(temp_re==0){memmove(x,xn1,sizeof(double)*4); return 0;}
    temp_re==1?(re=1):(re=re); 
    double xn2[4]={0};
    temp_re=cut_line_segment_from_line(xn1,l[2],xn2);
    //printf("panduan3 :%lf %lf %lf %lf\n",xn2[0],xn2[1],xn2[2],xn2[3]);

    if(temp_re==0){memmove(x,xn2,sizeof(double)*4); return 0;}
    temp_re==1?(re=1):(re=re);
    memmove(x,xn2,sizeof(double)*4);
    return re;
}




static inline int my_panduan2_is_in(double l[2][4],double* x)
{
    int re=0;
    double xn[4]={0};
    int temp_re=0;
    temp_re=cut_line_segment_from_line(x,l[0],xn);
    if(temp_re==0){return 0;}
    re=temp_re;
    double xn1[4]={0}; 
    temp_re=cut_line_segment_from_line(xn,l[1],xn1);
    if(temp_re==0){return 0;}
    temp_re==1?(re=1):(re=re); 
   
    return re;
}

//如果两个三角形共用一个顶点，那么直接判断一个三角形的边是否在角度范围内即可 

//
// 目前只支持凸多边形
//
static double** get_intersection_polygon_and_line(double** poly,int len,int *len1 ,double*l)
{        
    *len1=0;
    if(poly==NULL||len<=0)
    {
        return NULL;
    }
    double data[4]={0};
    data[0*2+0]=l[0];data[0*2+1]=l[1];
    data[1*2+0]=-l[1];data[1*2+1]=l[0];
    double data1[4]={0};
    data1[0*2+0]=l[0]; data1[0*2+1]=-l[1];
    data1[1*2+0]=l[1];data1[1*2+1]=l[0];

    double** temp_poly=(double**)malloc(sizeof(double*)*len);

    for(int i=0;i<len;i++)
    {
        double temp_p[2]={poly[i][0]-l[2],poly[i][1]-l[3]};
        temp_poly[i]=(double*)malloc(sizeof(double)*2);

        temp_poly[i][0]=data[0]*temp_p[0]+data[1]*temp_p[1];
        temp_poly[i][1]=data[2]*temp_p[0]+data[3]*temp_p[1];
        //printf("temp poly:%lf %lf\n",temp_poly[i][0],temp_poly[i][1]);
    }

    int bian_mark[2]={-1,-1};
    double p[2][2]={0},pp[2][2]={0};
    int k=0;
    for(int i=0;i<len;i++)
    {
        if(temp_poly[i][1]*temp_poly[(i+1)%len][1]<0)
        {
            bian_mark[k]=i;
            double lambda=temp_poly[(i+1)%len][1]/(temp_poly[(i+1)%len][1]-temp_poly[i][1]);
            p[k][0]=lambda*temp_poly[i][0]+(1-lambda)*temp_poly[(i+1)%len][0];
            p[k][1]=lambda*temp_poly[i][1]+(1-lambda)*temp_poly[(i+1)%len][1];

            pp[k][0]=data1[0]*p[k][0]+data1[1]*p[k][1];
            pp[k][1]=data1[2]*p[k][0]+data1[3]*p[k][1];
            pp[k][0]+=l[2];pp[k][1]+=l[3];
            k++;
        }
        if(k>=2){break;}
    }
    //printf("bian mark:%d %d k:%d\n",bian_mark[0],bian_mark[1],k);
    double ** re=NULL;
    if(k<=1)
    {
        if(temp_poly[0][1]>0)
        {
            re=(double**)malloc(sizeof(double*)*len);
            *len1=len;
            for(int i=0;i<len;i++)
            {
                re[i]=(double*)malloc(sizeof(double)*2);
                memmove(re[i],poly[i],sizeof(double)*2);
            }
        }    
    }
    else
    {
        re=(double**)malloc(sizeof(double*)*(len+1));
        //printf("temp poly:%lf \n",temp_poly[bian_mark[0]][1]);

        if(temp_poly[bian_mark[0]][1]<0)
        {
            *len1=0;
            re[*len1]=(double*)malloc(sizeof(double)*2);
            memmove(re[*len1],pp[0],sizeof(double)*2);
            *len1=(*len1+1);
            for(int i=0;i<len;i++)
            {
                re[*len1]=(double*)malloc(sizeof(double)*2);
                memmove(re[*len1],poly[(bian_mark[0]+i+1)%len],sizeof(double)*2);
                *len1=(*len1+1);
                if((bian_mark[0]+i+1)%len==bian_mark[1])
                {
                    re[*len1]=(double*)malloc(sizeof(double)*2);
                    memmove(re[*len1],pp[1],sizeof(double)*2);
                    *len1=(*len1+1);
                    break;
                }
            }
        }
        else
        {
            *len1=0;
            re[*len1]=(double*)malloc(sizeof(double)*2);
            memmove(re[*len1],poly[bian_mark[0]],sizeof(double)*2);
            *len1=(*len1+1);
            re[*len1]=(double*)malloc(sizeof(double)*2);
            memmove(re[*len1],pp[0],sizeof(double)*2);
            *len1=(*len1+1);

            re[*len1]=(double*)malloc(sizeof(double)*2);
            memmove(re[*len1],pp[1],sizeof(double)*2);
            *len1=(*len1+1);
            //printf("here len1:%d\n",*len1);
            for(int i=0;i<len;i++)
            {
                if((bian_mark[1]+i+1)%len==bian_mark[0])
                {
                    break;
                }
                re[*len1]=(double*)malloc(sizeof(double)*2);
                memmove(re[*len1],poly[(bian_mark[1]+i+1)%len],sizeof(double)*2);
                *len1=(*len1+1);
            }
        } 
    }

    for(int i=0;i<len;i++)
    {
        free(temp_poly[i]);
    }
    free(temp_poly);

    return re;
}

// 0 不相交
// 1 相交
//
static int is_two_triangle_intersecting(template_v* v1,template_v* v2,template_v*v3,
    template_v* v4,template_v* v5,template_v*v6)
{
    int num=0;
    template_v* vs[3]={NULL};
    int len=3;
    double** poly=(double**)malloc(sizeof(double*)*len);
    poly[0]=(double*)malloc(sizeof(double)*2);
    Greedy_Projection_Vert_Prop* gpvp1=(Greedy_Projection_Vert_Prop*)(v1->prop);
    memmove(poly[0],gpvp1->p,sizeof(double)*2);
    
    poly[1]=(double*)malloc(sizeof(double)*2);
    Greedy_Projection_Vert_Prop*gpvp2=(Greedy_Projection_Vert_Prop*)(v2->prop);
    memmove(poly[1],gpvp2->p,sizeof(double)*2);
    
    poly[2]=(double*)malloc(sizeof(double)*2);
    Greedy_Projection_Vert_Prop*gpvp3=(Greedy_Projection_Vert_Prop*)(v3->prop);
    memmove(poly[2],gpvp3->p,sizeof(double)*2);
    Greedy_Projection_Vert_Prop* gpvp4=(Greedy_Projection_Vert_Prop*)(v4->prop);
    Greedy_Projection_Vert_Prop* gpvp5=(Greedy_Projection_Vert_Prop*)(v5->prop);
    Greedy_Projection_Vert_Prop* gpvp6=(Greedy_Projection_Vert_Prop*)(v6->prop);
    if(v1==v4||v1==v5||v1==v6)
    {
        vs[num]=v1;num++;
    }
    if(v2==v4||v2==v5||v2==v6)
    {
        vs[num]=v2;num++;
    }
    if(v3==v4||v3==v5||v3==v6)
    {
        vs[num]=v3;num++;
    }
    int re=0;
    if(num==0)
    {
        //printf("num0\n");

        double l1[4]={gpvp5->p[0]-gpvp4->p[0],gpvp5->p[1]-gpvp4->p[1],gpvp4->p[0],gpvp4->p[1]};
        normalize(l1,2);
        double l2[4]={gpvp6->p[0]-gpvp5->p[0],gpvp6->p[1]-gpvp5->p[1],gpvp5->p[0],gpvp5->p[1]};
        normalize(l2,2);
        double l3[4]={gpvp4->p[0]-gpvp6->p[0],gpvp4->p[1]-gpvp6->p[1],gpvp6->p[0],gpvp6->p[1]};
        normalize(l3,2);
        double *l[3]={l1,l2,l3};
        for(int i=0;i<3;i++)
        {
            int len1=0;
            double**poly1= get_intersection_polygon_and_line(poly,len,&len1 ,l[i]);
            for(int i=0;i<len;i++)
            {
                free(poly[i]);
            }
            SAFE_FREE(poly);
            poly=poly1;len=len1;
            if(poly==NULL)
            {
                break;
            }
        } 
        poly==NULL?(re=0):(re=1);
        
    }
    else if(num==1)
    {
        //printf("num1\n");
        double l[2][4]={0};
        if(vs[0]==v1)
        {
            l[0][0]=gpvp2->p[0]-gpvp1->p[0];l[0][1]=gpvp2->p[1]-gpvp1->p[1];l[0][2]=gpvp1->p[0];l[0][3]=gpvp1->p[1];
            normalize(l[0],2);
            l[1][0]=gpvp1->p[0]-gpvp3->p[0];l[1][1]=gpvp1->p[1]-gpvp3->p[1];l[1][2]=gpvp1->p[0];l[1][3]=gpvp1->p[1];
            normalize(l[1],2);
        }
        else if(vs[0]==v2)
        {
            l[0][0]=gpvp3->p[0]-gpvp2->p[0];l[0][1]=gpvp3->p[1]-gpvp2->p[1];l[0][2]=gpvp2->p[0];l[0][3]=gpvp2->p[1];
            normalize(l[0],2);
            l[1][0]=gpvp2->p[0]-gpvp1->p[0];l[1][1]=gpvp2->p[1]-gpvp1->p[1];l[1][2]=gpvp2->p[0];l[1][3]=gpvp2->p[1];
            normalize(l[1],2);
        }
        else if(vs[0]==v3)
        {
            l[0][0]=gpvp1->p[0]-gpvp3->p[0];l[0][1]=gpvp1->p[1]-gpvp3->p[1];l[0][2]=gpvp3->p[0];l[0][3]=gpvp3->p[1];
            normalize(l[0],2);
            l[1][0]=gpvp3->p[0]-gpvp2->p[0];l[1][1]=gpvp3->p[1]-gpvp2->p[1];l[1][2]=gpvp3->p[0];l[1][3]=gpvp3->p[1];
            normalize(l[1],2);
        }
        double x[4]={0};
        if(vs[0]==v4)
        {
            x[0]=gpvp5->p[0];x[1]=gpvp5->p[1];x[2]=gpvp6->p[0];x[3]=gpvp6->p[1];
        }
        else if(vs[0]==v5)
        {
            x[0]=gpvp4->p[0];x[1]=gpvp4->p[1];x[2]=gpvp6->p[0];x[3]=gpvp6->p[1];
        }
        else if(vs[0]==v6)
        {
            x[0]=gpvp4->p[0];x[1]=gpvp4->p[1];x[2]=gpvp5->p[0];x[3]=gpvp5->p[1];
        }
        re=my_panduan2_is_in(l,x);
        re>0?(re=1):(re=re);
    }
    else
    {
        //printf("num2\n");

        re=0;
    }
    for(int i=0;i<len;i++)
    {
        free(poly[i]);
    }
    SAFE_FREE(poly);
    return re;
}

static inline int my_is_cell_legal1(template_c* c)
{
    template_v* v=NULL;
    template_hf*hf= (template_hf*)(c->halffaces->value);
    for(int i=0;i<c->vertices_size;i++)
    {
        if(c->vertices[i]!=hf->vertices[0]&&c->vertices[i]!=hf->vertices[1])
        {
            v=c->vertices[i];
            break;
        }
    }    
    Greedy_Projection_Vert_Prop* gpvp1=(Greedy_Projection_Vert_Prop*)(hf->vertices[0]->prop);
    Greedy_Projection_Vert_Prop* gpvp2=(Greedy_Projection_Vert_Prop*)(hf->vertices[1]->prop);
    Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
    
    double p1[2]={gpvp1->p[0]-gpvp->p[0],gpvp1->p[1]-gpvp->p[1]};
    double p2[2]={gpvp2->p[0]-gpvp->p[0],gpvp2->p[1]-gpvp->p[1]};

    if(p1[0]*p2[1]-p1[1]*p2[0]>0)
    {
        return 1;
    } 

    return 0;
}


static Node* touying_paichu(Int_RB_Tree* mark_hf, Node* node_hf,template_hf* hf,template_v* v)
{
    int jiance_id=94;
    if(v->id==jiance_id)
    {
        printf("begin once touying paichu %d %d\n",hf->vertices[0]->id,hf->vertices[1]->id);
    }
    if(mark_hf->find(mark_hf,hf->id)==NULL||node_hf==NULL)
    {
        if(v->id==jiance_id)
        {
            printf("end touying paichu1\n");
        }
        return node_hf;
    } 
    Node* re=NULL;
    Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop); 
    Greedy_Projection_Vert_Prop*gpvp1=(Greedy_Projection_Vert_Prop*)(hf->vertices[0]->prop);
    Greedy_Projection_Vert_Prop*gpvp2=(Greedy_Projection_Vert_Prop*)(hf->vertices[1]->prop);
    double p1[3]={gpvp1->p[0]-gpvp->p[0],gpvp1->p[1]-gpvp->p[1],0};
    double p2[3]={gpvp2->p[0]-gpvp->p[0],gpvp2->p[1]-gpvp->p[1],0};
    
    printf("gpvp:%lf %lf\n", gpvp->p[0],gpvp->p[1]);
    printf("gpvp1:%lf %lf\n",gpvp1->p[0],gpvp1->p[1]);
    printf("gpvp2:%lf %lf\n",gpvp2->p[0],gpvp2->p[1]);
 
    double l[3][4]={0};
    l[0][0]=gpvp1->p[0]-gpvp2->p[0];
    l[0][1]=gpvp1->p[1]-gpvp2->p[1];
    normalize(l[0],2);
    l[0][2]=gpvp1->p[0];l[0][3]=gpvp1->p[1];

    l[1][0]=p1[0];
    l[1][1]=p1[1];
    normalize(l[1],2);
    l[1][2]=gpvp1->p[0];l[1][3]=gpvp1->p[1];
        
    l[2][0]=-p2[0];
    l[2][1]=-p2[1];
    normalize(l[2],2);
    l[2][2]=gpvp2->p[0];l[2][3]=gpvp2->p[1]; 
    if(p1[0]*p2[1]-p1[1]*p2[0]<0)
    {
        l[0][0]=-l[0][0];l[0][1]=-l[0][1];
        l[1][0]=-l[1][0];l[1][1]=-l[1][1];
        l[2][0]=-l[2][0];l[2][1]=-l[2][1];
    }
    double xn[4]={0};
    for(Node* nit=node_hf;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf*hf1=(template_hf*)(nit->value);
        if(hf1==hf)
        {
            if(v->id==jiance_id)
            {
                printf("xiangdeng bu :%d %d\n",hf1->vertices[0]->id,hf1->vertices[1]->id);
            }
            re=node_overlying(re,hf1);
            continue;
        }
        
        gpvp1=(Greedy_Projection_Vert_Prop*)(hf1->vertices[0]->prop);
        gpvp2=(Greedy_Projection_Vert_Prop*)(hf1->vertices[1]->prop); 

        double* pp1=NULL;
        double *l1[2]={0};

        if(hf->vertices[0]==hf1->vertices[0] )
        {
            l1[0]=l[0];l1[1]=l[1];pp1=gpvp2->p;
        }
        else if(hf->vertices[0]==hf1->vertices[1])
        {
            l1[0]=l[0];l1[1]=l[1];pp1=gpvp1->p;
        }
        else if(hf->vertices[1]==hf1->vertices[0])
        {
            l1[0]=l[0];l1[1]=l[2];pp1=gpvp2->p;
        }
        else if(hf->vertices[1]==hf1->vertices[1])
        {
            l1[0]=l[0];l1[1]=l[2];pp1=gpvp1->p;
        }
        int mark=0; 
        if(pp1!=NULL)
        {
            my_panduan2_point_is_in(l1,pp1 )==1?(mark=2):(mark=0);
        }
        else
        {
            double x[4]={gpvp1->p[0],gpvp1->p[1],gpvp2->p[0],gpvp2->p[1] };

            if(v->id==jiance_id)
            {
                printf("x :%lf %lf %lf %lf\n",x[0],x[1],x[2],x[3]);
            }
            mark=my_panduan3_is_in(l,x);
            if(v->id==jiance_id)
            {
                printf("mark :%d v1:%d v2:%d\n",mark,hf1->vertices[0]->id,hf1->vertices[1]->id);
                printf("x :%lf %lf %lf %lf\n",x[0],x[1],x[2],x[3]);
            }
        }
        if(mark==2)
        {
            mark_hf->erase(mark_hf,hf1->id);
        }
        else if(mark==0)
        {           
            if(v->id==jiance_id)
            {
                printf("bu :%d %d\n",hf1->vertices[0]->id,hf1->vertices[1]->id);
            }

            re=node_overlying(re,hf1);
        } 
    }
    free_node(node_hf);
    if(v->id==jiance_id)
    {
        printf("end touying paichu\n");
    }
    return re;
}

static inline Node* chubu_paichu1(Mesh* m,double l[2][4],template_hf** hfs,Int_RB_Tree* mark_hf, Node* node_hf)
{
    Node* re=NULL;
    double x[4]={0};
    double *l1[2]={l[0],l[1]};
    for(Node* nit=node_hf;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_hf* hf=(template_hf*)(nit->value);
        Greedy_Projection_Vert_Prop* gpvp1=(Greedy_Projection_Vert_Prop*)(hf->vertices[0]->prop);
        Greedy_Projection_Vert_Prop* gpvp2=(Greedy_Projection_Vert_Prop*)(hf->vertices[1]->prop); 
        double* pp1=NULL;
        int mark=0;
        if(hf->vertices[0]==hfs[0]->vertices[1])
        {
            (my_panduan2_point_is_in(l1,gpvp2->p)==1)?(mark=2):(mark=0);
        }
        else if(hf->vertices[1]==hfs[0]->vertices[1])
        {
            (my_panduan2_point_is_in(l1,gpvp1->p)==1)?(mark=2):(mark=0);
        }
        else if(hf->vertices[0]==hfs[0]->vertices[0] )
        {
            (is_point_in_line_segment(l[0],gpvp2->p)==1)?(mark=2):(mark=0);
        }
        else if(hf->vertices[1]==hfs[0]->vertices[0])
        {
            (is_point_in_line_segment(l[0],gpvp1->p)==1)?(mark=2):(mark=0);
        }
        else if(hf->vertices[0]==hfs[1]->vertices[1])
        {
            (is_point_in_line_segment(l[1],gpvp2->p)==1)?(mark=2):(mark=0);
        } 
        else if(hf->vertices[1]==hfs[1]->vertices[1])
        {
            (is_point_in_line_segment(l[1],gpvp1->p)==1)?(mark=2):(mark=0);
        }
        else
        {
            x[0]=gpvp1->p[0];x[1]=gpvp1->p[1];
            x[2]=gpvp2->p[0];x[3]=gpvp2->p[1];
            mark= my_panduan2_is_in(l,x);
        }
        if(mark==2)
        {
            mark_hf->erase(mark_hf,hf->id);
        }
        else if(mark==0)
        {
            re=node_overlying(re,hf);
        }
    }
    return re;   
}

static inline Node* chubu_paichu(Mesh* m, template_v* v,Int_RB_Tree* mark_hf, Node* node_hf)
{
    Node* re=node_hf;
    template_hf* hfs[2]={0};
    double l[2][4]={0};
    for(auto vcit=m->vc_begin(m,*v);vcit!=m->vc_end(m,*v);vcit++)
    {
        if(!my_is_cell_legal1(quote(vcit)))
        {
            continue;
        }
        for(Node* nit=quote(vcit)->halffaces;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_hf* hf=(template_hf*)(nit->value);
            hf->vertices[1]==v?(hfs[0]=hf):(hf->vertices[0]==v?(hfs[1]=hf):(hfs[0]=hfs[0]) );
            
        }
        Greedy_Projection_Vert_Prop*gpvp1=(Greedy_Projection_Vert_Prop*)(hfs[0]->vertices[0]->prop);
        Greedy_Projection_Vert_Prop*gpvp2=(Greedy_Projection_Vert_Prop*)(hfs[0]->vertices[1]->prop); 
        l[0][0]=gpvp2->p[0]-gpvp1->p[0];l[0][1]=gpvp2->p[1]-gpvp1->p[1];l[0][2]=gpvp1->p[0];l[0][3]=gpvp1->p[1];
        normalize(l[0],2);
        gpvp1=(Greedy_Projection_Vert_Prop*)(hfs[1]->vertices[0]->prop);
        gpvp2=(Greedy_Projection_Vert_Prop*)(hfs[1]->vertices[1]->prop); 
        l[1][0]=gpvp2->p[0]-gpvp1->p[0];l[1][1]=gpvp2->p[1]-gpvp1->p[1];l[1][2]=gpvp1->p[0];l[1][3]=gpvp1->p[1];
        normalize(l[1],2);

        Node*temp_re= chubu_paichu1(m,l,hfs,mark_hf, re);
        free_node(re);
        re=temp_re;
    }
    return re;
}



//优先级排序
//创建的cell要和点法向同向

static inline void tinsert_one_vertex_to_mesh(Mesh* m,template_v* v,Int_RB_Tree* tree_c,Node** new_c)
{
    if(!m->vertex_is_boundary(m,*v))
    {
        return;
    }
    //template_c*c= vertex_is_in_which_cell(*new_c,v);
    template_c*c=vertex_is_in_which_cell_legal(*new_c,v);
    if(c!=NULL)
    {   if(v->cells!=NULL){return;}
        template_v* vs[3][3]={0};
        int i=0;
        for(Node* nit=c->halffaces;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_hf*hf=(template_hf*)(nit->value);
            vs[i][0]=hf->vertices[0];vs[i][1]=hf->vertices[1];vs[i][2]=v;
            i++;
        } 
        
        if(compute_score_triangle(vs[0][0],vs[0][1],vs[0][2])<=0||
            compute_score_triangle(vs[1][0],vs[1][1],vs[1][2])<=0||
            compute_score_triangle(vs[2][0],vs[2][1],vs[2][2])<=0 )
        {
            return;
        }
        *new_c=node_delete_value(*new_c,c);

        m->delete_cell(m,*c,true);
        template_c*c1= m->create_cellv(m,vs[0],3);
        c1!=NULL?(*new_c=node_overlying(*new_c,c1)):(*new_c=*new_c);
        if(c1!=NULL)
        {
            printf("create cell in:%d %d %d\n",c1->vertices[0]->id,c1->vertices[1]->id,c1->vertices[2]->id);
        }

        c1= m->create_cellv(m,vs[1],3);
        c1!=NULL?(*new_c=node_overlying(*new_c,c1)):(*new_c=*new_c);
        if(c1!=NULL)
        {
            printf("create cell in:%d %d %d\n",c1->vertices[0]->id,c1->vertices[1]->id,c1->vertices[2]->id);
        }
        c1=m->create_cellv(m,vs[2],3);
        c1!=NULL?(*new_c=node_overlying(*new_c,c1)):(*new_c=*new_c);
        if(c1!=NULL)
        {
            printf("create cell in:%d %d %d\n",c1->vertices[0]->id,c1->vertices[1]->id,c1->vertices[2]->id);
        }
        return ;
    }

    Node*node_hf=NULL,*node_all_hf=NULL;
    Int_RB_Tree* mark_hf=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(mark_hf);
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c1=(template_c*)(it.second);
        if(c1->vertices[0]==v||c1->vertices[1]==v||c1->vertices[2]==v)
        {
            continue;
        }
        for(auto chfit=m->chf_begin(m,*c1);chfit!=m->chf_end(m,*c1);chfit++)
        {
            if(m->face_is_boundary(m,*(quote(chfit)->face)))
            {
                mark_hf->insert(mark_hf,quote(chfit)->id,quote(chfit));
                Greedy_Projection_Vert_Prop* gpvp1=(Greedy_Projection_Vert_Prop*)(quote(chfit)->vertices[0]->prop);
                Greedy_Projection_Vert_Prop* gpvp2=(Greedy_Projection_Vert_Prop*)(quote(chfit)->vertices[1]->prop); 
                
                if(gpvp1->active==1&&gpvp2->active==1)
                {
                    node_hf=node_overlying(node_hf,quote(chfit));
                } 
            }
        }
    }
    for(Node* nit=*new_c;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c1=(template_c*)(nit->value);
        for(auto chfit=m->chf_begin(m,*c1);chfit!=m->chf_end(m,*c1);chfit++)
        {
            if(c1->vertices[0]==v||c1->vertices[1]==v||c1->vertices[2]==v)
            {
                continue;
            }

            if(m->face_is_boundary(m,*(quote(chfit)->face)))
            {
                node_hf=node_overlying(node_hf,quote(chfit));
                mark_hf->insert(mark_hf,quote(chfit)->id,quote(chfit));
            }
        }
    }

    printf("node hf size:%d\n",node_size(node_hf));
    node_hf=chubu_paichu(m,v,mark_hf,node_hf);

    printf("node hf size:%d\n",node_size(node_hf));

    // for(Node* nit=node_hf;nit!=NULL;nit=(Node*)(nit->Next))
    // {
    //     template_hf* hf=(template_hf*)(nit->value);
    //     printf("hf v1:%d v2:%d\n",hf->vertices[0]->id,hf->vertices[1]->id);
    // }

    for(auto it=mark_hf->begin(mark_hf);it.it!=NULL;it++)
    {
        node_all_hf=node_overlying(node_all_hf,it.second);
    }

    for(Node* nit=node_all_hf;nit!=NULL;nit=(Node*)(nit->Next))
    {
        node_hf=touying_paichu(mark_hf, node_hf,(template_hf*)(nit->value),v);
    }
    printf("node hf size:%d\n",node_size(node_hf));

    for(Node*nit=node_hf;nit!=NULL;nit=(Node*)(nit->Next) )
    {
        template_hf*hf=(template_hf*)(nit->value);
        Greedy_Projection_Vert_Prop* gpvp1=(Greedy_Projection_Vert_Prop*)(hf->vertices[1]->prop);
        Greedy_Projection_Vert_Prop* gpvp2=(Greedy_Projection_Vert_Prop*)(hf->vertices[0]->prop);
        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        double p1[3]={gpvp1->p[0]-gpvp->p[0],gpvp1->p[1]-gpvp->p[1]};
        double p2[3]={gpvp2->p[0]-gpvp->p[0],gpvp2->p[1]-gpvp->p[1]};
        //static int is_legal_triangle_angle(double *p1,double* p2,double* p3 )
        //is_legal_triangle_angle(v->point,hf->vertices[1]->point,hf->vertices[0]->point)
        if(p1[0]*p2[1]-p1[1]*p2[0]>0&&compute_score_triangle(hf->vertices[1],hf->vertices[0],v)>0 )
        {
            template_v*vs[3]={hf->vertices[1],hf->vertices[0],v};
            template_c*c1= m->create_cellv(m,vs,3);
            if(c1!=NULL)
            {
                printf("create cell out:%d %d %d\n",c1->vertices[0]->id,c1->vertices[1]->id,c1->vertices[2]->id);
            }
            c1!=NULL?(*new_c=node_overlying(*new_c,c1)):(*new_c=*new_c);
        }
    }
    free_node(node_all_hf);
    int_rb_tree_free(mark_hf);
    free_node(node_hf);
}

static void fixed_cells_intersecting_tree_cn(Mesh*m,Int_RB_Tree* tree_c )
{   
    printf("begin fix cells\n");
    if(tree_c->size<=0){return ;}
    int len=tree_c->size;
    int *mark1=(int*)malloc(sizeof(int)*len);
    memset(mark1,0,sizeof(int)*len);
    int** mark=(int**)malloc(sizeof(int*)*len);
    template_c** cells=(template_c**)malloc(sizeof(template_c*)*len);
    
    int temp_i=0;
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        cells[temp_i]=(template_c*)(it.second);
        temp_i++;
    }

    for(int i=0;i<len;i++)
    {
        mark[i]=(int*)malloc(sizeof(int)*len);
        memset(mark[i],0,sizeof(int)*len);
    }
   
    int flag=1;
    while(flag)
    {
        temp_i=-1;
        for(int i=0;i<len;i++)
        {   
            if(cells[i]==NULL||mark1[i]==1)
            {continue;}
            if(m->cell_is_boundary(m,*(cells[i])))
            {
                for(int j=0;j<len;j++)
                {
                    if(i==j||cells[j]==NULL||mark1[j]==1||mark[i][j]==1){continue;}
                    if(is_two_triangle_intersecting(cells[j]->vertices[0],cells[j]->vertices[1],cells[j]->vertices[2],
                        cells[i]->vertices[0],cells[i]->vertices[1],cells[i]->vertices[2]))
                    {
                        temp_i=i;break;
                    }
                    else
                    {
                        mark[i][j]=1;mark[j][i]=1;
                    }
                }
                if(temp_i!=-1){break;} 
                else{mark1[i]=1;}
            }   
        }
        if(temp_i==-1){flag=0;}
        else
        {
            tree_c->erase(tree_c,cells[temp_i]->id);
            m->delete_cell(m,*(cells[temp_i]),true);
            cells[temp_i]=NULL;
        }
    }
    free(mark1);
    free(cells);
    for(int i=0;i<len;i++)
    {
        free(mark[i]);
    }
    free(mark);
    printf("end fix cells\n");
}

static Int_RB_Tree* get_legal_cells_tree_c(Mesh* m,Int_RB_Tree* tree_c,int is_fix)
{
    Int_RB_Tree* re=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(re);
    int mark=0;
    if(is_fix)
    {
        Node* temp_node=NULL;
        for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
        {
            template_c* c=(template_c*)(it.second);
            Greedy_Projection_Vert_Prop* gpvp1=(Greedy_Projection_Vert_Prop*)(c->vertices[0]->prop);
            Greedy_Projection_Vert_Prop* gpvp2=(Greedy_Projection_Vert_Prop*)(c->vertices[1]->prop);
            Greedy_Projection_Vert_Prop* gpvp3=(Greedy_Projection_Vert_Prop*)(c->vertices[2]->prop);
            
            if(compute_score_triangle(c->vertices[0],c->vertices[1],c->vertices[2])<=0)
            {
                temp_node=node_overlying(temp_node,c);
                continue;
            }
            if(gpvp1->active==1&&gpvp2->active==1&&gpvp3->active==1)
            {
                if(!my_is_cell_legal1(c))
                {
                    temp_node=node_overlying(temp_node,c);
                    continue;
                }
                else
                {
                    c->prop=&mark;
                }
            }
        }

        for(Node*nit=temp_node;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_c* c=(template_c*)(nit->value);
            tree_c->erase(tree_c,c->id);
            m->delete_cell(m,*c,true);
        }
        free_node(temp_node);
    }
    
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        if(c->prop!=NULL)
        {
            c->prop=NULL;
            re->insert(re,c->id,c);
            continue;
        }
        if(my_is_cell_legal1(c))
        {
            re->insert(re,c->id,c);
        }
    }
    if(is_fix)
    {
        fixed_cells_intersecting_tree_cn(m,re);
    }
    return re;
}


static inline int my_panduan_f_legal(template_f*f)
{
    return (f->halffaces[0]->cell->prop!=NULL&&f->halffaces[1]->cell->prop!=NULL);
}
static void tiaozheng_new_cellsn(Mesh* m,Node*new_c,Int_RB_Tree* tree_c)
{
    
    Int_RB_Tree * tree_f=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_f);
    int mark=0;
    for(Node* nit=new_c;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        c->prop=&mark;
    }
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        c->prop=&mark;
    }

    for(Node* nit=new_c;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        for(auto chfit=m->chf_begin(m,*c);chfit!=m->chf_end(m,*c);chfit++)
        {
            if(!m->face_is_boundary(m,*(quote(chfit)->face)))
            {
                if(my_panduan_f_legal(quote(chfit)->face))
                {
                    tree_f->insert(tree_f,quote(chfit)->face->id,quote(chfit)->face);
                }
            }
        }
    }

    for(Node* nit=new_c;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_c* c=(template_c*)(nit->value);
        c->prop=NULL;
    }
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        c->prop=NULL;
    }
    Int_RB_Tree * temp_mark=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(temp_mark);
    int temp_num=0;
    int flag=1;
    while(flag)
    {
        printf("once temp_num:%d***************************************************\n",temp_num);
        Int_RB_Tree * tree_f1=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
        int_rb_tree_init(tree_f1);
        flag=0;
        for(auto it=tree_f->begin(tree_f);it.it!=NULL;it++)
        {
            //这里判断
            template_f*f=m->get_facep(m,it.first);
            if(f==NULL||m->face_is_boundary(m,*f))
            {
                continue;
            }
            if(temp_mark->find(temp_mark,f->id)!=NULL)
            {
                tree_f1->insert(tree_f1,f->id,f);
                continue;
            }
            template_f*f1= touying_tiaozheng_edge(m, f);
            //还要确保f的两个cell都是tree_c1成员
            //确保libcell操作更加纯净
            //确保f 没有被删除
            //
            if(f1!=NULL)
            {
                tree_f1->insert(tree_f1,f1->id,f1);flag=1;
            }
            else 
            {
                f=m->get_facep(m,it.first);
                if(f!=NULL)
                {
                    temp_mark->insert(temp_mark,f->id,f);
                    tree_f1->insert(tree_f1,f->id,f);
                }
            }
        }  
        
        temp_num++;
        int_rb_tree_free(tree_f);
        tree_f=tree_f1;

        if(temp_num>40)
        {
            break;
        }
    }
    int_rb_tree_free(temp_mark);
    int_rb_tree_free(tree_f);
}

static int my_cmp(void* a,void*b)
{
    template_v* va= (template_v*)a,* vb= (template_v*)b;
    Greedy_Projection_Vert_Prop*gpvp1=(Greedy_Projection_Vert_Prop*)(va->prop);
    Greedy_Projection_Vert_Prop*gpvp2=(Greedy_Projection_Vert_Prop*)(vb->prop);

    return  gpvp1->theta>gpvp2->theta;

}
// 
void tanlan_touying1(template_v* v,double* normal, Node* node_v,Mesh* m,double r,int is_fix );

void tanlan_touying2(template_v* v,double* normal, Node* node_v,Mesh* m,double r)
{
    Node* temp_node_v=NULL;
    for(Node* nit1=node_v;nit1!=NULL;nit1=(Node*)(nit1->Next))
    {
        template_v* v=(template_v*)(nit1->value);
        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        if(inner_product(normal,gpvp->normal,3)>0)
        {
            temp_node_v=node_overlying(temp_node_v,v); 
        }
    }

    for(Node* nit=temp_node_v;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        gpvp->theta=normal[0]*gpvp->normal[0]+ normal[1]*gpvp->normal[1]+normal[2]*gpvp->normal[2]; 
    }
    printf("once \n");
    Node* node_v1=node_quick_sort(temp_node_v,my_cmp);
    // for(Node* nit=node_v1;nit!=NULL;nit=(Node*)(nit->Next))
    // {
    //     template_v* v=(template_v*)(nit->value);
    //     Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
    //     printf(" v id:%d theta : %lf\n",v->id,gpvp->theta);
    // }

    tanlan_touying1( v, normal, node_v1,m,r,1);
   
    for(Node* nit=temp_node_v;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        gpvp->theta=0; 
    }
    free_node(node_v1);
    free_node(temp_node_v);
}

static inline int int_rb_tree_cell_is_empty(Mesh* m,Int_RB_Tree* tree_c)
{
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {   

        template_c* c1=(template_c*)(it.second);
        
        for(auto chfit=m->chf_begin(m,*c1);chfit!=m->chf_end(m,*c1);chfit++)
        {
            if(m->face_is_boundary(m,*(quote(chfit)->face)))
            {
                Greedy_Projection_Vert_Prop* gpvp1=(Greedy_Projection_Vert_Prop*)(quote(chfit)->vertices[0]->prop);
                Greedy_Projection_Vert_Prop* gpvp2=(Greedy_Projection_Vert_Prop*)(quote(chfit)->vertices[1]->prop);
                if(gpvp1->active==1&&gpvp2->active==1)
                {
                    return 0;
                } 
            }
        }
    }
    return 1;
}

// static int is_two_triangle_intersecting(template_v* v1,template_v* v2,template_v*v3,
//     template_v* v4,template_v* v5,template_v*v6)

static inline int int_rb_tree_cell_is_intersecting(Int_RB_Tree* tree_c,
    template_v* v1,template_v* v2,template_v* v3)
{
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        if(is_two_triangle_intersecting(v1,v2,v3,
            c->vertices[0],c->vertices[1],c->vertices[2]))
        {
            return 1;
        }
    } 

    return 0;
}


// //此参数的normal不是v的normal，而应该是node_v的平均normal
// tree_c中排除不合适的cell,也就是法向相反的cell
// //首先排除node_v中法向相反的点
// // 首先排除node_v 中落在初始tree_c中的点
//  
void tanlan_touying1(template_v* v,double* normal, Node* node_v,Mesh* m,double r,int is_fix)
{ 
    double theta = acos(normal[2]);
    double  n1[3]={normal[1],-normal[0],0};
    double*matd = get_matrix_from_rotation(theta,n1);
    
    Int_RB_Tree * tree_c=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree_c);
    for(Node*nit=node_v;nit!=NULL;nit=(Node*)(nit->Next) )
    {
        template_v* v1=(template_v*)(nit->value);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v1->prop);
        gpvp->is_init=1;memset(gpvp->p,0,sizeof(double)*3);gpvp->theta=0;gpvp->is_legal=0;gpvp->active=1; 
        for(auto vcit=m->vc_begin(m,*v1);vcit!=m->vc_end(m,*v1);vcit++ )
        {
            if(tree_c->insert(tree_c,quote(vcit)->id,quote(vcit)))
            {
                for(int k=0;k<3;k++)
                {
                    gpvp=(Greedy_Projection_Vert_Prop*)(quote(vcit)->vertices[k]->prop);
                    gpvp->is_init=1;memset(gpvp->p,0,sizeof(double)*3);gpvp->theta=0;gpvp->is_legal=0;  
                }  
            }
        }    
    }
    for(Node* nit=node_v;nit!=NULL;nit=(Node*)(nit->Next))
    {
        touying_one_vertex(matd,v,(template_v*)(nit->value)); 
    }
    
    Node* tool_node=NULL;
    for(auto it=tree_c->begin(tree_c);it.it!=NULL;it++)
    {
        template_c* c=(template_c*)(it.second);
        for(int i=0;i<3;i++)
        {
            Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(c->vertices[i]->prop);
            if(gpvp->is_init==1)
            {
                tool_node=node_overlying(tool_node,c->vertices[i]);
                touying_one_vertex(matd,v,c->vertices[i]); 
            }   
        } 
    }

    Int_RB_Tree * tree_c1=get_legal_cells_tree_c(m,tree_c,is_fix);


    Node* new_c=NULL;
    Node* node_v1=NULL;
    for(Node*nit=node_v;nit!=NULL;nit=(Node*)(nit->Next) )
    {
        template_v* v1=(template_v*)(nit->value);
        if(m->vertex_is_boundary(m,*v1)&&vertex_is_in_which_cell1(tree_c1,v1)==NULL)
        { 
            node_v1=node_overlying(node_v1,v1);
        }
    }
    
    Node* nnit=node_v1;
    //
    //这里要修改
    //要改为第一个node_hf为NULL
    //且新建的cell不和tree_c1相交
    //
    if(int_rb_tree_cell_is_empty(m,tree_c1))
    {
        if(node_size(node_v1)<3)
        {
            goto myend;
        }
        template_v* vs[3]={0};
        for(int i=0;i<3;i++)
        {
            vs[i]=(template_v*)(nnit->value);
            nnit=(Node*)(nnit->Next);  
        }

        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(vs[0]->prop);

        Greedy_Projection_Vert_Prop* gpvp1=(Greedy_Projection_Vert_Prop*)(vs[1]->prop);
        Greedy_Projection_Vert_Prop* gpvp2=(Greedy_Projection_Vert_Prop*)(vs[2]->prop);
        
        double p1[3]={gpvp1->p[0]-gpvp->p[0],gpvp1->p[1]-gpvp->p[1]};
        double p2[3]={gpvp2->p[0]-gpvp->p[0],gpvp2->p[1]-gpvp->p[1]};
        template_c* c1=NULL;
        if(p1[0]*p2[1]-p1[1]*p2[0]<0)
        {
            template_v* temp_v=vs[0];
            vs[0]=vs[1];
            vs[1]=temp_v;
        }

        if(int_rb_tree_cell_is_intersecting(tree_c1,
            vs[0],vs[1],vs[2]))
        {
            goto myend;
        }
        c1=m->create_cellv(m,vs,3);
        if(c1!=NULL)
        {
            new_c=node_overlying(new_c,c1);
            printf("create init cell:%d %d %d\n",c1->vertices[0]->id,c1->vertices[1]->id,c1->vertices[2]->id);
        }

    } 
    for(;nnit!=NULL;nnit=(Node*)(nnit->Next) )
    {
        template_v* v1=(template_v*)(nnit->value);
        printf("begin t insert v id:%d\n",v1->id);
        tinsert_one_vertex_to_mesh(m,v1,tree_c1,&new_c);
    }

    tiaozheng_new_cellsn(m,new_c,tree_c1);
 

// static int touying_tiaozheng_edge(Mesh* m, template_f* f)

// 
myend:

    for(Node*nit=node_v;nit!=NULL;nit=(Node*)(nit->Next) )
    {
        template_v* v1=(template_v*)(nit->value);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v1->prop);
        gpvp->is_init=1;memset(gpvp->p,0,sizeof(double)*3);gpvp->theta=0;gpvp->is_legal=0;gpvp->active=0;
    } 
    for(Node*nit=tool_node;nit!=NULL;nit=(Node*)(nit->Next) )
    {
        template_v* v1=(template_v*)(nit->value);
        Greedy_Projection_Vert_Prop*gpvp=(Greedy_Projection_Vert_Prop*)(v1->prop);
        gpvp->is_init=1;memset(gpvp->p,0,sizeof(double)*3);gpvp->theta=0;gpvp->is_legal=0;gpvp->active=0;
    }
 

    free_node(tool_node);
    free_node(new_c);
    free_node(node_v1);
    int_rb_tree_free(tree_c);
    int_rb_tree_free(tree_c1);
    SAFE_FREE(matd);

}
typedef struct My_Visual_ChongJian{
    Node* opnit; 
    Viewer_Something* vsp;
    Viewer_Something* vsf;
    Viewer_Something* vse;
    Viewer_Opengl_Interpreter* voi;

}My_Visual_ChongJian;

static inline void my_visual_chongjian_init(My_Visual_ChongJian* mvcj )
{
    mvcj->opnit=NULL; 
    mvcj->vsp=NULL;
    mvcj->vsf=NULL;
    mvcj->vse=NULL;
    mvcj->voi=NULL;

}
static void my_visual_chongjian_key_callback(Viewer_Intera*vi)
{
    Interactor_GlobalInfo* g_info=vi->g_info;
    My_Visual_ChongJian* mvcj=(My_Visual_ChongJian*)(vi->representation);

;
    if(g_info->key_action==VIEWER_PRESS&&g_info->key==VIEWER_KEY_F)
    {
        if(mvcj->opnit==NULL||mvcj->opnit->Next==NULL)
        {
            return;
        }
        mvcj->opnit=(Node*)(mvcj->opnit->Next);
        Points_Open_Set* pos=(Points_Open_Set*)(mvcj->opnit->value);
        Viewer_Points* vp=(Viewer_Points*)(mvcj->vsp->evolution);
        vp->Data_rows=node_size(pos->n);
        SAFE_FREE(vp->Data);
        vp->Data=(float*)malloc(sizeof(float)*3*vp->Data_rows);
        int i=0;
        for(Node* nit=pos->n;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_v* v=(template_v*)(nit->value);
            vp->Data[i*3+0]=v->point[0];vp->Data[i*3+1]=v->point[1];vp->Data[i*3+2]=v->point[2];
            i++;
        }
        vp->set_color(vp,1.0,1.0,0.0,1.0);
        
        Viewer_Edges* ve=(Viewer_Edges*)(mvcj->vse->evolution);
        double normal[3]={0};
        for(Node* nit1=pos->n;nit1!=NULL;nit1=(Node*)(nit1->Next))
        {
            template_v* v=(template_v*)(nit1->value);
            Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
            normal[0]+=gpvp->normal[0];normal[1]+=gpvp->normal[1];normal[2]+=gpvp->normal[2];
        }
        normalize(normal,3);

        ve->Data[0]=pos->v->point[0];ve->Data[1]=pos->v->point[1];ve->Data[2]=pos->v->point[2];
        ve->Data[3]=pos->v->point[0]+normal[0]*0.5;ve->Data[4]=pos->v->point[1]+normal[1]*0.5;ve->Data[5]=pos->v->point[2]+normal[2]*0.5;
    
        ve->Data_index[0]=0;ve->Data_index[1]=1; 

        mvcj->voi->update_data(mvcj->voi);
    }
    else if(g_info->key_action==VIEWER_PRESS&&g_info->key==VIEWER_KEY_R)
    {
        if(mvcj->opnit==NULL||mvcj->opnit->Prev==NULL)
        {
            return;
        }
        mvcj->opnit=(Node*)(mvcj->opnit->Prev);
        Points_Open_Set* pos=(Points_Open_Set*)(mvcj->opnit->value);
        Viewer_Points* vp=(Viewer_Points*)(mvcj->vsp->evolution);
        vp->Data_rows=node_size(pos->n);
        SAFE_FREE(vp->Data);
        vp->Data=(float*)malloc(sizeof(float)*3*vp->Data_rows);
        int i=0;
        for(Node* nit=pos->n;nit!=NULL;nit=(Node*)(nit->Next))
        {
            template_v* v=(template_v*)(nit->value);
            vp->Data[i*3+0]=v->point[0];vp->Data[i*3+1]=v->point[1];vp->Data[i*3+2]=v->point[2];
            i++;
        }
        vp->set_color(vp,1.0,1.0,0.0,1.0);
        
        Viewer_Edges* ve=(Viewer_Edges*)(mvcj->vse->evolution);
        double normal[3]={0};
        for(Node* nit1=pos->n;nit1!=NULL;nit1=(Node*)(nit1->Next))
        {
            template_v* v=(template_v*)(nit1->value);
            Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
            normal[0]+=gpvp->normal[0];normal[1]+=gpvp->normal[1];normal[2]+=gpvp->normal[2];
        }
        normalize(normal,3);

        ve->Data[0]=pos->v->point[0];ve->Data[1]=pos->v->point[1];ve->Data[2]=pos->v->point[2];
        ve->Data[3]=pos->v->point[0]+normal[0]*0.5;ve->Data[4]=pos->v->point[1]+normal[1]*0.5;ve->Data[5]=pos->v->point[2]+normal[2]*0.5;
    
        ve->Data_index[0]=0;ve->Data_index[1]=1; 

        mvcj->voi->update_data(mvcj->voi);
    }
}
void test_touying_chongjian()
{
    Mesh mesh;
    Mesh_init(&mesh);
    _ReadOff_(&mesh,"lawen1.off",3);
    for(auto cit=mesh.c_begin(&mesh);cit!=mesh.c_end(&mesh);cit++)
    {
        quote(cit)->prop=compute_cell_normal(quote(cit));
    }
    Node* node_v=NULL;
    Mesh mesh1;
    Mesh_init(&mesh1);
    mesh1.dimension=2;mesh1.simplex=1; mesh1.manifold_require=1;
     
    for(auto vit=mesh.v_begin(&mesh);vit!=mesh.v_end(&mesh);vit++)
    {
        // double* normal=(double*)malloc(sizeof(double)*3);
        // memset(normal,0,sizeof(double)*3);
        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)malloc(sizeof(Greedy_Projection_Vert_Prop));
        greedy_projection_vert_prop_init(gpvp);
        for(auto vcit=mesh.vc_begin(&mesh,*vit);vcit!=mesh.vc_end(&mesh,*vit);vcit++)
        {
            gpvp->normal[0]+=((double*)(quote(vcit)->prop))[0];
            gpvp->normal[1]+=((double*)(quote(vcit)->prop))[1];
            gpvp->normal[2]+=((double*)(quote(vcit)->prop))[2];
        }
        normalize(gpvp->normal,3);
        template_v* v= mesh1.create_vertexv(&mesh1,quote(vit)->point,3);
        v->prop=gpvp;
        node_v=node_overlying(node_v,v);
        // if(normalize(normal,3))
        // {
        //     template_v* v= mesh1.create_vertexv(&mesh1,quote(vit)->point,3);
        //     v->prop=normal;
        //     node_v=node_overlying(node_v,v);
        //     //quote(vit)->prop=normal;
        // }
    }

    Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    int_rb_tree_init(tree);
    printf("begin compute open sets\n");

    template_v**values=(template_v**)malloc(sizeof(template_v*)*node_size(node_v));
    int i=0;

    for(Node* nit=node_v;nit!=NULL;nit=(Node*)(nit->Next))
    {
        values[i]=(template_v*)(nit->value);
        i++;
    }
    KD_Node* kdtree= create_kd_tree(values,i,0);
       
    Node*poses= compute_points_open_sets(node_v,tree,0.20,kdtree);
    printf("end compute open sets\n");
    int chongjian_i=0;Node* poses_nit=poses;
    for(Node* nit=poses;nit!=NULL;nit=(Node*)(nit->Next))
    {
        printf("chongjian_i :%d *******************************************\n",chongjian_i);

        poses_nit=nit;
        Points_Open_Set* pos=(Points_Open_Set*)(nit->value);
        double normal[3]={0};
        for(Node* nit1=pos->n;nit1!=NULL;nit1=(Node*)(nit1->Next))
        {
            template_v* v=(template_v*)(nit1->value);
            Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
            normal[0]+=gpvp->normal[0];normal[1]+=gpvp->normal[1];normal[2]+=gpvp->normal[2];
        }
        normalize(normal,3);
        printf("end temp v\n");
        // tanlan_touying1(pos->v,normal, temp_node_v,&mesh1,1.5);
        tanlan_touying2(pos->v,normal, pos->n,&mesh1,1.5);

        if(chongjian_i>=21960)
        {
            break;
        }
        chongjian_i++;
    } 
    template_v* vv_z=mesh1.get_vertexp(&mesh1,69933);
    Node* nearest_s= kd_tree_find_nearest_sphere(vv_z,1.5*0.2,kdtree);
    printf("begin nearest_s********************************begin nearest_s\n");

    double normaln[3]={0};
    for(Node* nit1=nearest_s;nit1!=NULL;nit1=(Node*)(nit1->Next))
    {
        template_v* v=(template_v*)(nit1->value);
        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        normaln[0]+=gpvp->normal[0];normaln[1]+=gpvp->normal[1];normaln[2]+=gpvp->normal[2];
    }
    normalize(normaln,3);
    tanlan_touying2(vv_z,normaln, nearest_s,&mesh1,1.5);
    free_node(nearest_s);

    vv_z=mesh1.get_vertexp(&mesh1,69920);
    nearest_s= kd_tree_find_nearest_sphere(vv_z,1.5*0.2,kdtree);
    printf("begin nearest_s********************************begin nearest_s\n");
    memset(normaln,0,sizeof(double)*3);
    for(Node* nit1=nearest_s;nit1!=NULL;nit1=(Node*)(nit1->Next))
    {
        template_v* v=(template_v*)(nit1->value);
        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        normaln[0]+=gpvp->normal[0];normaln[1]+=gpvp->normal[1];normaln[2]+=gpvp->normal[2];
    }
    normalize(normaln,3);
    tanlan_touying2(vv_z,normaln, nearest_s,&mesh1,1.5);
    free_node(nearest_s);
    vv_z=mesh1.get_vertexp(&mesh1,69121);
    nearest_s= kd_tree_find_nearest_sphere(vv_z,1.5*0.2,kdtree);
    printf("begin nearest_s********************************begin nearest_s\n");
    memset(normaln,0,sizeof(double)*3);
    for(Node* nit1=nearest_s;nit1!=NULL;nit1=(Node*)(nit1->Next))
    {
        template_v* v=(template_v*)(nit1->value);
        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        normaln[0]+=gpvp->normal[0];normaln[1]+=gpvp->normal[1];normaln[2]+=gpvp->normal[2];
    }
    normalize(normaln,3);
    tanlan_touying2(vv_z,normaln, nearest_s,&mesh1,1.5);
    free_node(nearest_s);

    vv_z=mesh1.get_vertexp(&mesh1,69212);
    nearest_s= kd_tree_find_nearest_sphere(vv_z,1.5*0.2,kdtree);
    printf("begin nearest_s********************************begin nearest_s\n");
    memset(normaln,0,sizeof(double)*3);
    for(Node* nit1=nearest_s;nit1!=NULL;nit1=(Node*)(nit1->Next))
    {
        template_v* v=(template_v*)(nit1->value);
        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        normaln[0]+=gpvp->normal[0];normaln[1]+=gpvp->normal[1];normaln[2]+=gpvp->normal[2];
    }
    normalize(normaln,3);
    tanlan_touying2(vv_z,normaln, nearest_s,&mesh1,1.5);
    free_node(nearest_s);

    vv_z=mesh1.get_vertexp(&mesh1,69117);
    nearest_s= kd_tree_find_nearest_sphere(vv_z,1.5*0.2,kdtree);
    printf("begin nearest_s********************************begin nearest_s\n");
    memset(normaln,0,sizeof(double)*3);
    for(Node* nit1=nearest_s;nit1!=NULL;nit1=(Node*)(nit1->Next))
    {
        template_v* v=(template_v*)(nit1->value);
        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        normaln[0]+=gpvp->normal[0];normaln[1]+=gpvp->normal[1];normaln[2]+=gpvp->normal[2];
    }
    normalize(normaln,3);
    tanlan_touying2(vv_z,normaln, nearest_s,&mesh1,1.5);
    free_node(nearest_s);
    vv_z=mesh1.get_vertexp(&mesh1,69099);
    nearest_s= kd_tree_find_nearest_sphere(vv_z,1.5*0.2,kdtree);
    printf("begin nearest_s********************************begin nearest_s\n");
    memset(normaln,0,sizeof(double)*3);
    for(Node* nit1=nearest_s;nit1!=NULL;nit1=(Node*)(nit1->Next))
    {
        template_v* v=(template_v*)(nit1->value);
        Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
        normaln[0]+=gpvp->normal[0];normaln[1]+=gpvp->normal[1];normaln[2]+=gpvp->normal[2];
    }
    normalize(normaln,3);
    tanlan_touying2(vv_z,normaln, nearest_s,&mesh1,1.5);
    free_node(nearest_s);
    printf("node pose size:%d\n",node_size(poses));
    // template_v* ver=mesh1.get_vertexp(&mesh1,12);
    // printf("%lf %lf %lf\n",ver->point[0],ver->point[1],ver->point[2]); 
// Points name_id:0 something_id:0 marked_id: 12
// piont :-0.504755 -25.722561 -0.797140


    // Mesh mesh2;
    // Mesh_init(&mesh2);


    Viewer_World_Manager vwm;
    viewer_world_manager_init(&vwm);
    Viewer_World *vw=vwm.create_world(&vwm,NULL);
    Viewer_Opengl_Interpreter voi;

    Viewer_Arcroll*va=add_default_somethings(vw);


    Node* n=vw->create_something(vw,"Points");
    Viewer_Something* vs=(Viewer_Something*)(n->value);
    //vs->disappear=1;
    Viewer_Points* vp=(Viewer_Points*)(vs->evolution);
    
    vp->Data_rows=node_size(node_v);
    vp->Data=(float*)malloc(sizeof(float)*3*vp->Data_rows);
    i=0;
    for(Node* nit=node_reverse(node_v);nit!=NULL;nit=(Node*)(nit->Prev))
    {
        template_v* v=(template_v*)(nit->value);
        vp->Data[i*3+0]=v->point[0];vp->Data[i*3+1]=v->point[1];vp->Data[i*3+2]=v->point[2];
        i++;
    }
    vp->pointsize=6.0;
    vp->set_color(vp,1.0,0.0,0.0,1.0);
    free_node(n); 
    //***************
    n=vw->create_something(vw,"Faces");    
    vs=(Viewer_Something*)(n->value);
   // ms->disappear=1;
    Viewer_Faces* vf=(Viewer_Faces*)(vs->evolution);
    

    vf->Data_rows=mesh1.num_v(&mesh1);
    vf->Data_index_rows=mesh1.num_c(&mesh1);
    vf->color_rows=vf->Data_rows;
    vf->normal_rows=vf->Data_index_rows;
    vf->set_color(vf,1.0,1.0,1.0,1.0);

    get_data_from_2dim_cell(&mesh1,&(vf->Data),&(vf->Data_index));
        //vf->is_reversal_normal=1;
 
    free_node(n);
    //******************
//*****************
    // n=vw->create_something(vw,"Edges");
    // vs=(Viewer_Something*)(n->value);
    // Viewer_Edges* ve=(Viewer_Edges*)(vs->evolution);
    // ve->Data_rows=node_size(node_v1)*2;
    // ve->Data_index_rows=ve->Data_rows;
    // ve->Data=(float*)malloc(sizeof(float)*3*ve->Data_rows);
    // ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows);
    // i=0;
    // for(Node* nit=node_v1;nit!=NULL;nit=(Node*)(nit->Next))
    // {
    //     template_v* v=(template_v*)(nit->value);
    //     ve->Data[i*3+0]=v->point[0];ve->Data[i*3+1]=v->point[1];ve->Data[i*3+2]=v->point[2];
    //     double* normal=(double*)(v->prop);
    //     i++;
    //     ve->Data[i*3+0]=v->point[0]+0.1*normal[0];ve->Data[i*3+1]=v->point[1]+0.1*normal[1];ve->Data[i*3+2]=v->point[2]+0.1*normal[2];
    //     i++;
    // }
    // for(i=0;i<ve->Data_index_rows;i++)
    // {
    //     ve->Data_index[i*2+0]=i*2+0 ; ve->Data_index[i*2+1]= i*2+1;
    // }
    // ve->color_rows=ve->Data_index_rows;
    // ve->set_color(ve,1.0,1.0,0.0,1.0);
    // free_node(n);  
    //**********************
    My_Visual_ChongJian* mvcj=(My_Visual_ChongJian*)malloc(sizeof(My_Visual_ChongJian));
    my_visual_chongjian_init(mvcj);
    Points_Open_Set* pos=(Points_Open_Set*)(poses->value);
    n=vw->create_something(vw,"Edges");
    vs=(Viewer_Something*)(n->value);
    Viewer_Edges* ve=(Viewer_Edges*)(vs->evolution);
    ve->Data_rows=2;
    ve->Data_index_rows=1;
    ve->Data=(float*)malloc(sizeof(float)*3*ve->Data_rows);
    ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows);

    double normal[3]={0};
    for(Node* nit1=pos->n;nit1!=NULL;nit1=(Node*)(nit1->Next))
    {
            template_v* v=(template_v*)(nit1->value);
            Greedy_Projection_Vert_Prop* gpvp=(Greedy_Projection_Vert_Prop*)(v->prop);
            normal[0]+=gpvp->normal[0];normal[1]+=gpvp->normal[1];normal[2]+=gpvp->normal[2];
    }
    normalize(normal,3);

    ve->Data[0]=pos->v->point[0];ve->Data[1]=pos->v->point[1];ve->Data[2]=pos->v->point[2];
    ve->Data[3]=pos->v->point[0]+normal[0]*0.5;ve->Data[4]=pos->v->point[1]+normal[1]*0.5;ve->Data[5]=pos->v->point[2]+normal[2]*0.5;
    
    ve->Data_index[0]=0;ve->Data_index[1]=1;
    ve->color_rows=ve->Data_index_rows;
    ve->set_color(ve,0,1.0,1.0,1.0);


    free_node(n);
   
// Viewer_Something* viewer_tools_create_3d_text(wchar_t *str,float* xyz,LB_Matrix*mat,Viewer_World*vw );
// void viewer_tools_append_3d_text(wchar_t *str,float* xyz,LB_Matrix*mat,Viewer_World*vw,Viewer_Something* vsi );


    mvcj->vse=vs;

    ThreeD_Text_Intera* tdti=(ThreeD_Text_Intera*)malloc(sizeof(ThreeD_Text_Intera));
    threed_text_intera_init(tdti);
    wchar_t temp_str[40]={0};
    //threed_tex 
    n=vw->create_something(vw,"Points");
    vs=(Viewer_Something*)(n->value);
    vp=(Viewer_Points*)(vs->evolution);
    vp->pointsize=9.0;
    vp->Data_rows=node_size(  ((Points_Open_Set*)(poses_nit->value))->n  );
    vp->Data=(float*)malloc(sizeof(float)*3*vp->Data_rows);
    i=0;
    for(Node* nit=((Points_Open_Set*)(poses_nit->value))->n;nit!=NULL;nit=(Node*)(nit->Next))
    {

        template_v* v=(template_v*)(nit->value);
        
        vp->Data[i*3+0]=v->point[0];vp->Data[i*3+1]=v->point[1];vp->Data[i*3+2]=v->point[2];
        swprintf(temp_str,40,L"%d",v->id);
        viewer_tools_append_3d_textn(temp_str,&(vp->Data[i*3+0]),vp->mat,vw,tdti);

        i++;
    }

    vp->set_color(vp,1.0,1.0,0.0,1.0);
    free_node(n);  

    Viewer_Something*vsi= viewer_tools_create_3d_textn(vw,tdti );


   
    mvcj->opnit=poses;
    mvcj->vsp=vs;
    mvcj->voi=&voi;
      
    n=vw->create_something(vw,"Intera");
    vs=(Viewer_Something*)(n->value);
    Viewer_Intera*vi=(Viewer_Intera*)(vs->evolution);
    vi->representation=(void*)mvcj;
    // vi->cursor_position_callback=viewer_Arcroll_cursor_position_callback;
    // vi->scroll_callback=viewer_Arcroll_scroll_callback;
    // vi->mouse_button_callback=viewer_Arcroll_mouse_button_callback;
    vi->key_callback=my_visual_chongjian_key_callback; 
    // //vi->drop_callback= viewer_Arcroll_drop_callback;
    // //vi->animation=viewer_Arcroll_animation;
    // //vi->char_callback=viewer_Arcroll_char_callback;
    free_node(n); 


    viewer_opengl_interpreter_initn(&voi,&vwm);

    voi.interpreter(&voi);
    viewer_opengl_interpreter_free(&voi);
    viewer_world_manager_free(&vwm);


    free_kdnode(kdtree);

    free_node(poses);

    int_rb_tree_free(tree);
    Mesh_free(&mesh);
}

void test_sampling()
{
    Mesh mesh;
    Mesh_init(&mesh);

    _ReadOff_(&mesh,"lawen.off",3);
    mesh.external_cell_init_(&mesh);

    Node* node_v=create_surface_uniform_sampling_points(&mesh,30000,0);

    for(Node* nit=node_v;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        template_c*c=(template_c*)(v->prop);
        if(c!=NULL)
        {
            Cell_Prop* cp=(Cell_Prop*)(c->prop);
            if(cp==NULL)
            {
                cp=(Cell_Prop*)malloc(sizeof(Cell_Prop));
                cell_prop_init(cp);
                cp->normal=compute_cell_normal(c);
                c->prop=cp;
            }
            cp->node=node_overlying(cp->node,v);
        }
    }
    template_v** values=(template_v**)malloc(sizeof(template_v*)*node_size(node_v));
    int i=0;
    for(Node* nit=node_v;nit!=NULL;nit=(Node*)(nit->Next))
    {
        values[i]=(template_v*)(nit->value); 
        i++;
    } 
    //printf("begin create kd_ tree\n"); 
    KD_Node* kn=create_kd_tree(values,i,0);
    free(values);

    printf("node _v:%d\n",node_size(node_v));
    double hou=0.0;
    Node* node_v1=NULL;

    for(Node* nit=node_v;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        template_c*c=(template_c*)(v->prop);
        if(c!=NULL)
        {
            Cell_Prop* cp=(Cell_Prop*)(c->prop);
            Vertex*vn=(Vertex*)malloc(sizeof(Vertex));
            Vertex_init_(vn);
            vn->point_size=3;
            vn->point=(double*)malloc(sizeof(double)*3);
            vn->point[0]=v->point[0]+hou*cp->normal[0] ;
            vn->point[1]=v->point[1]+hou*cp->normal[1];
            vn->point[2]=v->point[2]+hou*cp->normal[2] ;

            // KD_Node*kdn= kd_tree_find_nearest(vn,kn);
            // if(distance_of_two_points(vn->point,kdn->value->point,3 )<0.9*hou)
            // {
            //     free_Vertex(vn);
            // }
            // else
            {
                node_v1=node_overlying(node_v1,vn);
                vn->prop=v;
            }

        }
    }
    values=(template_v**)malloc(sizeof(template_v*)*node_size(node_v1));
    i=0;
    for(Node* nit=node_v1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        values[i]=(template_v*)(nit->value); 
        i++;
    } 
    KD_Node*kdtree1= create_kd_tree(values,i,0); 
    free(values);
    my_points_cloud_compute_normalr(node_v1,kdtree1,0.4);
    
    free_kdnode(kdtree1);
   

    // Mesh mesh1;
    // Mesh_init(&mesh1);
    // mesh1.dimension=2;mesh1.simplex=0;mesh1.manifold_require=1;
    // double pp[3]={0};
    // Int_RB_Tree* tree=(Int_RB_Tree*)malloc(sizeof(Int_RB_Tree));
    // int_rb_tree_init(tree);
    // Node* node_hfs=NULL;
    // for(auto fit=mesh.f_begin(&mesh);fit!=mesh.f_end(&mesh);fit++)
    // {
    //     if(mesh.face_is_boundary(&mesh,*fit))
    //     {
    //         if(quote(fit)->halffaces[0]->cell!=NULL)
    //         {
    //             node_hfs= node_overlying(node_hfs,quote(fit)->halffaces[0]);
    //         }
    //         else if(quote(fit)->halffaces[1]->cell!=NULL)
    //         {
    //             node_hfs=node_overlying(node_hfs,quote(fit)->halffaces[1]);
    //         }
    //         template_v* v1=(template_v*)(tree->find(tree,quote(fit)->vertices[0]->id));
    //         template_v* v2=(template_v*)(tree->find(tree,quote(fit)->vertices[1]->id));
    //         //double p[3]={0};
    //         if(v1==NULL)
    //         {
    //             v1=compute_hou_point_from_vertex(&mesh,quote(fit)->vertices[0],hou,&mesh1);
    //             tree->insert(tree,quote(fit)->vertices[0]->id,v1);
    //         }
    //         if(v2==NULL)
    //         {
    //             v2=compute_hou_point_from_vertex(&mesh,quote(fit)->vertices[1],hou,&mesh1);
    //             tree->insert(tree,quote(fit)->vertices[1]->id,v2);
    //         }
    //         template_v* vs[2]={v1,v2};   
    //         mesh1.create_facev(&mesh1,vs,2);
    //     }
    // }  

    // for(auto vit=mesh1.v_begin(&mesh1);vit!=mesh1.v_end(&mesh1);vit++)
    // {
    //     KD_Node*kdn= kd_tree_find_nearest(quote(vit),kn);
    //     if(distance_of_two_points(quote(vit)->point,kdn->value->point,3 )<0.95*hou)
    //     {
    //         template_v* vs[2]={NULL};
    //         template_f* fs[2]={NULL};
    //         int j=0,k=0;
    //         for(auto vfit=mesh1.vf_begin(&mesh1,*vit);vfit!=mesh1.vf_end(&mesh1,*vit);vfit++)
    //         {
    //             fs[k]=quote(vfit);k++;
    //             if(quote(vfit)->vertices[0]!=quote(vit))
    //             {
    //                 vs[j]=quote(vfit)->vertices[0]; j++;
    //             }
    //             else if(quote(vfit)->vertices[1]!=quote(vit) )
    //             {
    //                 vs[j]=quote(vfit)->vertices[1]; j++;
    //             }
    //         } 
    //         mesh1.delete_face(&mesh1,*(fs[0]),1);
    //         mesh1.delete_face(&mesh1,*(fs[1]),1);
    //         printf("delete once\n");
    //         if(j==2)
    //         {
    //             mesh1.create_facev(&mesh1,vs,2);
    //         }
    //     }
    // }  
    

    // for(Node* nit=mesh.external_cell.halffaces;nit!=NULL;nit=(Node*)(nit->Next))
    // {
    //     template_hf*hf=(template_hf*)(nit->value);
    //     template_v* v1=hf->vertices[1],*v2=hf->vertices[0],*v3=NULL;
    //     if(nit->Next==NULL)
    //     {
    //         v3=((template_hf*)(mesh.external_cell.halffaces->value))->vertices[0];
    //     }
    //     else
    //     {
    //         v3=((template_hf*)(((Node*)(nit->Next))->value))->vertices[0]; 
    //     }
    //     double normal[3]={0};
    //     for(auto vcit=mesh.vc_begin(&mesh,*v2);vcit!=mesh.vc_end(&mesh,*v2);vcit++)
    //     {
    //         Cell_Prop* cp=(Cell_Prop*)(quote(vcit)->prop);
    //         if(cp!=NULL)
    //         {
    //             normal[0]+=cp->normal[0];
    //             normal[1]+=cp->normal[1];
    //             normal[2]+=cp->normal[2];
    //         }
    //     } 
    //     if(normalize(normal,3))
    //     {
    //         double p1[3]={v1->point[0]-v2->point[0],v1->point[1]-v2->point[1],v1->point[2]-v2->point[2]};
    //         double p2[3]={v3->point[0]-v2->point[0],v3->point[1]-v2->point[1],v3->point[2]-v2->point[2]};
    //         double inn1=inner_product(p1,normal,3);
    //         double inn2=inner_product(p2,normal,3);
    //         p1[0]-=normal[0]*inn1;p1[1]-=normal[1]*inn1;p1[2]-=normal[2]*inn1;
    //         p2[0]-=normal[0]*inn2;p2[1]-=normal[1]*inn2;p2[2]-=normal[2]*inn2;
    //         normalize(p1,3); normalize(p2,3);
    //         double p4[3]={(p1[0]+p2[0])/2.0,(p1[1]+p2[1])/2.0,(p1[2]+p2[2])/2.0};
    //         double* p3=out_product(p1,p2);
    //         if(inner_product(normal,p3,3)>0)
    //         {
    //             p4[0]=-p4[0];p4[1]=-p4[1];p4[2]=-p4[2];
    //         }
    //         normalize(p4,3);
    //         SAFE_FREE(p3);
    //         double p5[3]={v2->point[0]+hou*normal[0],v2->point[1]+hou*normal[1],
    //             v2->point[2]+hou*normal[2]};
    //         template_v* vn=  mesh1.create_vertexv(&mesh1,p5,3);
    //         // Vertex*vn=(Vertex*)malloc(sizeof(Vertex));
    //         // Vertex_init_(vn);
    //         // vn->point_size=3;
    //         // vn->point=(double*)malloc(sizeof(double)*3);
    //         // vn->point[0]=v2->point[0]+hou*normal[0];
    //         // vn->point[1]=v2->point[1]+hou*normal[1];
    //         // vn->point[2]=v2->point[2]+hou*normal[2];
    //         double* normal1=(double*)malloc(sizeof(double)*3);
    //         normal1[0]=normal[0];normal1[1]=normal[1];normal1[2]=normal[2];
    //         vn->prop=normal1;
    //         //node_v1=node_overlying(node_v1,vn);
    //         tree->insert(tree,v2->id,vn);
    //     }
    // }
    // mesh1.dimension=2;mesh1.simplex=0;mesh1.manifold_require=1;
    
    // for(Node* nit=mesh.external_cell.halffaces;nit!=NULL;nit=(Node*)(nit->Next))
    // {
    //     template_hf*hf=(template_hf*)(nit->value);
    //     template_v* v1=hf->vertices[1],*v2=hf->vertices[0];
    //     template_v* vn1=(template_v*)(tree->find(tree,v1->id));
    //     template_v* vn2=(template_v*)(tree->find(tree,v2->id));
    //     template_v* vs[2]={vn1,vn2};
    //     if(vn1!=NULL&&vn2!=NULL)
    //     {
    //         mesh1.create_facev(&mesh1,vs,2);
    //     }
    // }
    // int_rb_tree_free(tree);

    // for(Node* nit=mesh.external_cell.halffaces;nit!=NULL;nit=(Node*)(nit->Next))
    // {
    //     template_hf*hf=(template_hf*)(nit->value);
    //     template_v* v1=hf->vertices[1],*v2=hf->vertices[0],*v3=NULL;
    //     //printf("%d %d\n",hf->vertices[0]->id,hf->vertices[1]->id );

    //     if(nit->Next==NULL)
    //     {
    //         v3=((template_hf*)(mesh.external_cell.halffaces->value))->vertices[0];
    //     }
    //     else
    //     {
    //         v3=((template_hf*)(((Node*)(nit->Next))->value))->vertices[0]; 
    //     }
    //     double normal[3]={0};
    //     for(auto vcit=mesh.vc_begin(&mesh,*v2);vcit!=mesh.vc_end(&mesh,*v2);vcit++)
    //     {
    //         Cell_Prop* cp=(Cell_Prop*)(quote(vcit)->prop);
    //         if(cp!=NULL)
    //         {
    //             normal[0]+=cp->normal[0];
    //             normal[1]+=cp->normal[1];
    //             normal[2]+=cp->normal[2];
    //         }
    //     } 
    //     if(normalize(normal,3))
    //     {
    //         double p1[3]={v1->point[0]-v2->point[0],v1->point[1]-v2->point[1],v1->point[2]-v2->point[2]};
    //         double p2[3]={v3->point[0]-v2->point[0],v3->point[1]-v2->point[1],v3->point[2]-v2->point[2]};
    //         double inn1=inner_product(p1,normal,3);
    //         double inn2=inner_product(p2,normal,3);
    //         p1[0]-=normal[0]*inn1;p1[1]-=normal[1]*inn1;p1[2]-=normal[2]*inn1;
    //         p2[0]-=normal[0]*inn2;p2[1]-=normal[1]*inn2;p2[2]-=normal[2]*inn2;
    //         normalize(p1,3); normalize(p2,3);

    //         double p4[3]={(p1[0]+p2[0])/2.0,(p1[1]+p2[1])/2.0,(p1[2]+p2[2])/2.0};

    //         double* p3=out_product(p1,p2);

    //         if(inner_product(normal,p3,3)<0)
    //         {
    //             p4[0]=-p4[0];p4[1]=-p4[1];p4[2]=-p4[2];
    //         }
    //         normalize(p4,3);

    //         SAFE_FREE(p3);
    //         for( i=0;i<5;i++)
    //         {
    //             Vertex*vn=(Vertex*)malloc(sizeof(Vertex));
    //             Vertex_init_(vn);
    //             vn->point_size=3;
    //             vn->point=(double*)malloc(sizeof(double)*3);
    //             vn->point[0]=v2->point[0]+hou*normal[0]*i/5.0;
    //             vn->point[1]=v2->point[1]+hou*normal[1]*i/5.0;
    //             vn->point[2]=v2->point[2]+hou*normal[2]*i/5.0;
    //             double* normal1=(double*)malloc(sizeof(double)*3);
    //             normal1[0]=p4[0];normal1[1]=p4[1];normal1[2]=p4[2];
    //             vn->prop=normal1;
    //             node_v1=node_overlying(node_v1,vn);
    //         }
    //     }
    // }
   
    //node_v1=node_v;

    double* box=mesh_compute_bounding_box_from_node(node_v1);


    double centroid[3]={(box[0]+box[3])/2.0,(box[1]+box[4])/2.0,(box[2]+box[5])/2.0};
    printf("centroid : %lf %lf %lf\n",centroid[0],centroid[1],centroid[2]);
    double max=box[3]-box[0];

    max=((box[4]-box[1])>max?((box[4]-box[1])):(max));
    max=((box[5]-box[2])>max?((box[5]-box[2])):(max));   
    printf("max:%lf\n",max);
    for(Node* nit=node_v1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v*v=(template_v*)(nit->value);
        v->point[0]-=centroid[0];
        v->point[1]-=centroid[1];
        v->point[2]-=centroid[2];
        v->point[0]*= (1.9/max);
        v->point[1]*= (1.9/max);
        v->point[2]*= (1.9/max);
    }
    // for(auto vit=mesh1.v_begin(&mesh1);vit!=mesh1.v_end(&mesh1);vit++)
    // {
    //     quote(vit)->point[0]-=centroid[0];
    //     quote(vit)->point[1]-=centroid[1];
    //     quote(vit)->point[2]-=centroid[2];
    //     quote(vit)->point[0]*= (1.95/max);
    //     quote(vit)->point[1]*= (1.95/max);
    //     quote(vit)->point[2]*= (1.95/max); 

    // }

    free(box); 
   
    Mesh* m1=create_poisson_reconstruct_mesh(node_v1,310,500);

    // Poisson_Grid*pg=(Poisson_Grid*)malloc(sizeof(Poisson_Grid));
    // poisson_grid_init(pg);

    // poisson_grid_init_with_depth(pg,7);

    // fill_normal_coeff(pg,node_v1);

    // poisson_test_sparse_matrix(pg);
    
    // double iso=compute_average_poisson_iso(pg);
    // printf("begin hrerere\n");
    // // Node* node_tris=marching_cubes_generate_triangles(0.2,100,my_fun2,NULL);
    // // Node* node_tris=marching_cubes_generate_triangles(iso,100,end_poisson_space_fun,pg);
    // Node*node_tris1=NULL;

    // //Node****node_tris=marching_cubes_generate_trianglesn(iso,100,end_poisson_space_fun,pg,&node_tris1);

    // Node****node_tris=my_marching_cubes_generate_trianglesn(node_v1,iso,150,end_poisson_space_fun,pg,&node_tris1);

    // Mesh* m1= create_mesh_from_marching_cubes(node_tris,node_tris1,150);

    
    //double r=compute_visual_avarage_distance(node_v1);

    //adjust_normal_of_node_vertices(node_v1,r,kdtree1);


    Viewer_World_Manager vwm;
    viewer_world_manager_init(&vwm);
    Viewer_World *vw=vwm.create_world(&vwm,NULL);

    add_default_somethings(vw);  
    //test_show_triangles(vw,node_tris1);

    test_show_mesh_cells(vw,m1);

    Node* n=vw->create_something(vw,"Points");
    Viewer_Something* vs=(Viewer_Something*)(n->value);
    Viewer_Points* vp=(Viewer_Points*)(vs->evolution);
    
    vp->Data_rows=node_size(node_v);
    vp->Data=(float*)malloc(sizeof(float)*3*vp->Data_rows);
    i=0;
    for(Node* nit=node_v;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        vp->Data[i*3+0]=v->point[0];vp->Data[i*3+1]=v->point[1];vp->Data[i*3+2]=v->point[2];
        i++;
    }

    vp->set_color(vp,1.0,0.0,0.0,1.0);
    free_node(n); 
    
    n=vw->create_something(vw,"Points");
    vs=(Viewer_Something*)(n->value);
    vp=(Viewer_Points*)(vs->evolution);
    
    vp->Data_rows=node_size(node_v1);
    vp->Data=(float*)malloc(sizeof(float)*3*vp->Data_rows);
    i=0;
    for(Node* nit=node_v1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        vp->Data[i*3+0]=v->point[0];vp->Data[i*3+1]=v->point[1];vp->Data[i*3+2]=v->point[2];
        i++;
    }

    vp->set_color(vp,0.0,1.0,1.0,1.0);
    free_node(n);

    // n=vw->create_something(vw,"Edges");
    // vs=(Viewer_Something*)(n->value);
    // Viewer_Edges* ve=(Viewer_Edges*)(vs->evolution);
    // ve->Data_rows=mesh1.num_v(&mesh1);
    // ve->Data_index_rows=mesh1.num_f(&mesh1);

    // get_lines_data_from_2dim_cell(&mesh1 ,&(ve->Data) ,&(ve->Data_index));

    // ve->color_rows=ve->Data_index_rows;
    // ve->set_color(ve,1.0,0.0,0.0,1.0);
    // free_node(n); 

    n=vw->create_something(vw,"Edges");
    vs=(Viewer_Something*)(n->value);
    Viewer_Edges* ve=(Viewer_Edges*)(vs->evolution);

    ve->Data_rows=node_size(node_v1)*2;
    ve->Data_index_rows=node_size(node_v1);

    ve->Data=(float*)malloc(sizeof(float)*3*ve->Data_rows);
    ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows);

    i=0;
    for(Node* nit=node_v1;nit!=NULL;nit=(Node*)(nit->Next))
    {
        template_v* v=(template_v*)(nit->value);
        ve->Data[i*3+0]=v->point[0];ve->Data[i*3+1]=v->point[1];ve->Data[i*3+2]=v->point[2];
        double* normal=(double*)(v->prop);
        i++;
        ve->Data[i*3+0]=v->point[0]+0.1*normal[0];ve->Data[i*3+1]=v->point[1]+0.1*normal[1];ve->Data[i*3+2]=v->point[2]+0.1*normal[2];

        i++;
    }
    for(i=0;i<ve->Data_index_rows;i++)
    {
        ve->Data_index[i*2+0]=i*2+0 ; ve->Data_index[i*2+1]= i*2+1;
    }
    ve->color_rows=ve->Data_index_rows;
    ve->set_color(ve,1.0,1.0,0.0,1.0);
    free_node(n); 

    //******************
    // n=vw->create_something(vw,"Edges");
    // vs=(Viewer_Something*)(n->value);
    // Viewer_Edges* ve=(Viewer_Edges*)(vs->evolution);

    // ve->Data_rows=node_size(node_v)*2;
    // ve->Data_index_rows=node_size(node_v);

    // ve->Data=(float*)malloc(sizeof(float)*3*ve->Data_rows);
    // ve->Data_index=(unsigned int*)malloc(sizeof(unsigned int)*2*ve->Data_index_rows);

    // i=0;
    // for(Node* nit=node_v;nit!=NULL;nit=(Node*)(nit->Next))
    // {
    //     template_v* v=(template_v*)(nit->value);
    //     ve->Data[i*3+0]=v->point[0];ve->Data[i*3+1]=v->point[1];ve->Data[i*3+2]=v->point[2];
    //     template_c*c=(template_c*)(v->prop);
    //     Cell_Prop*cp=(Cell_Prop*)(c->prop);
    //     i++;
    //     ve->Data[i*3+0]=v->point[0]+cp->normal[0];ve->Data[i*3+1]=v->point[1]+cp->normal[1];ve->Data[i*3+2]=v->point[2]+cp->normal[2];

    //     i++;
    // }
    // for(i=0;i<ve->Data_index_rows;i++)
    // {
    //     ve->Data_index[i*2+0]=i*2+0 ; ve->Data_index[i*2+1]= i*2+1;
    // }
    // ve->color_rows=ve->Data_index_rows;
    // ve->set_color(ve,1.0,1.0,0.0,1.0);
    // free_node(n); 
 
    Viewer_Opengl_Interpreter voi;
    viewer_opengl_interpreter_initn(&voi,&vwm);
    voi.interpreter(&voi);
    viewer_opengl_interpreter_free(&voi);
    viewer_world_manager_free(&vwm);

    Mesh_free(&mesh);
}

int main()
{
    printf("%lf\n",atan2(-1,-1));
    test_touying_chongjian();
    // double p1[2]={0,0},p2[2]={1,0},p3[2]={0.9,1.9},p4[2]={0,1};
    // double iscon= is_concyclic_four_vertices(p2,p3,p1,p4);
    // printf("is con :%lf\n",iscon);

    double p1[3]={0,0,0},p2[3]={1,0,0},p3[3]={0,1,0};
    double q1[3]={0.5,1,0},q2[3]={0,3.0,0},q3[3]={-0.5,-0.5,0};
    double* poly[3]={p1,p2,p3};
    int len=3;
    int len1=0;
    double l[4]={q2[0]-q1[0],q2[1]-q1[1],q1[0],q1[1]};
    normalize(l,2);

    double**poly1= get_intersection_polygon_and_line(poly,len,&len1 ,l);


    double l1[4]={q3[0]-q2[0],q3[1]-q2[1],q2[0],q2[1]};
    normalize(l1,2);
    int len2=0;
    double**poly2= get_intersection_polygon_and_line(poly1,len1,&len2 ,l1);

    double l2[4]={q1[0]-q3[0],q1[1]-q3[1],q3[0],q3[1]};
    normalize(l2,2);
    int len3=0;
    double**poly3= get_intersection_polygon_and_line(poly2,len2,&len3 ,l2);

    printf("len1:%d\n",len3);
    for(int i=0;i<len3;i++)
    {
        printf("%lf %lf\n",poly3[i][0],poly3[i][1]);
    }



    //test_sampling();
    //test_poisson();
    printf("end\n");

	return 0;
}